module util_cmfd
    implicit none
!    real(8), pointer, dimension(:,:) :: hDif, hTot, hRmv, hNuF
!    real(8), pointer, dimension(:,:,:) :: hmSct
    integer, private, parameter :: fdhat = 777, fdtilde = 877
    integer, private, parameter :: fj=91, fdif=92, fphiac=93, fcmfd=94, fMTCMFD=95, fMTCMFDOFF=96, fRL1CMFD=97, fRL2CMFD=98, ft=501
    integer, private, parameter :: fcj=1091, fcp=1093
    integer, private, parameter :: fblc=260
    integer, private :: fcmfdf(100)
    real(8), pointer, dimension(:, :) :: phi
    real(8), private, pointer, dimension(:) :: psi, psip
    real(8) :: l2ref, sumref
    real(8), private :: tally_log
    
    real(8), private, pointer, dimension(:,:) :: gHF
    real(8), private, pointer, dimension(:,:,:) :: gHSCT
    real(8), private :: keff
    
    real(8), private, pointer, dimension(:,:,:) :: gHF3
    
    
!    integer, private :: stcnt
!    real(8), private, pointer, dimension(:) :: sumpsi, sqsumpsi
    integer, private :: ndim
    
    real(8), private, pointer, dimension(:,:) :: fsigt, fsigtb
    
contains

subroutine prepare_cmfd
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "CMFD_LARSEN.FI"
    integer :: g, i
    character(128) :: fn
    

    
!    open(1234, file='dbg.txt', status='unknown')

!    allocate(hDif(totm_cmfd, ng))
!    allocate(hTot(totm_cmfd, ng))
!    allocate(hRmv(totm_cmfd, ng))
!    allocate(hNuF(totm_cmfd, ng))
!    allocate(hmSct(ng,ng, totm_cmfd))
!    allocate(phiftemp(totmesh, ng))
    ndim=totm_cmfd
    allocate(phi(totm_cmfd, ng))
    allocate(psi(totm_cmfd))
    allocate(psip(totm_cmfd))
    
    
    ! gcmfd3
    allocate(gHF(0:totm_cmfd+1,ng))
    allocate(gHSCT(ng,ng,0:totm_cmfd+1))
    allocate(gHF3(3,0:totm_cmfd+1,ng))
    
!    allocate(sumpsi(totm_cmfd))
!    allocate(sqsumpsi(totm_cmfd))
!    stcnt = 0
    
    phi = 1.
    
    if (PDFREF) then
        sumref = sum(refc)
        l2ref = 0
        do g=1, ng
            do i=1, totm_cmfd
                l2ref = l2ref + refc(i, g)**2
            end do
        end do
        l2ref = sqrt(l2ref)
    end if 
    
#ifdef _PRTDTLD
    open(fdtilde, file='dtildeo.txt', status='unknown')
#endif
!    open(fdhat, file='dhato.txt', status='unknown')
!    open(fdhat+1, file='dhate.txt', status='unknown')
!    open(fdhat+2, file='dhat.txt', status='unknown')
        
!    open(fj, file='NetJ.txt', status='unknown')
!    open(fdif, file='Diff.txt', status='unknown')
!    open(fphiac, file='CMFDfluxac.txt', status='unknown')
!    open(fcj, file='CJ.txt', status='unknown')
!    open(fcp, file='CF.txt', status='unknown')
    open(fcmfd, file='CMFD.txt', status='unknown')
    open(ft, file="CMFDACC.txt", status='unknown')
    
    open(fblc, file="nbalance.txt", status='unknown')
    
    
!    do g=1, ng
!        write(fn, '(a,I<MAX(0,0)>, a)'), "CMFD_FLUX", g, "G.txt"
!        fcmfdf(g)=3000+g
!        open(fcmfdf(g), file=fn, status='unknown')
!    enddo
!    if (ng .eq. 2) open(fcmfd2, file='CMFD2.txt', status='unknown')
!    open(fMTCMFD, file='MT_CMFD.txt', status='unknown')
!    open(fMTCMFDOFF, file='MT_CMFDOFF.txt', status='unknown')
!    open(fRL1CMFD, file='MT_R1CMFD.txt', status='unknown')
!    open(fRL2CMFD, file='MT_R2CMFD.txt', status='unknown')
    
!    open(1001, file="CHK_DIF.txt", status='unknown')
    
    
!#ifdef CMFDLARSEN
!    dtildeL(0) = 1._8/3/Xtot(complist(1),1)/delx_2(0)
!    do i=1, totm_cmfd-1
!        dtildeL(i) =(h_cmfd(i)+h_cmfd(i+1))/(Xtot(complist(i),1)*h_cmfd(i) + Xtot(complist(i+1),1)*h_cmfd(i+1))/3/delx_2(i) 
!    enddo
!    dtildeL(i) = 1._8/3/Xtot(complist(i),1)/delx_2(i)
!#endif
    if (cmfdtype .eq. 7 .or. cmfdtype .eq. 9) then
        allocate(fsigt(totm_cmfd,ng))
        allocate(fsigtb(totm_cmfd,ng))
        fsigt=0
        fsigtb=0
        call eval_sigtbar()
    endif
end subroutine

subroutine eval_sigtbar
    include "pspec.FI"
    include "xsec.FI"
    integer :: r, i, ii, is, ie, g
    real(8) :: hj, xj
    r = totmesh/totm_cmfd
    ! evaluate fsigt and fsigtb
    do g=1, ng; do i=1, totm_cmfd-1
        ! left
        hj=h_cmfd(i)
        xj=hac_cmfd(i-1)
        is=r*(i-1)+1
        ie=r*i
        do ii=is, ie
            fsigt(i,g)=fsigt(i,g)+1._8/(2*hj)*( (hac(ii)-xj)**2 - (hac(ii-1)-xj)**2 )*Xtot(complist(ii),g)
        enddo
        ! right
        hj=h_cmfd(i+1)
        xj=hac_cmfd(i+1)
        is=r*i+1
        ie=r*(i+1)
        do ii=is, ie
            fsigt(i,g)=fsigt(i,g)-1._8/(2*hj)*( (hac(ii)-xj)**2 - (hac(ii-1)-xj)**2 )*Xtot(complist(ii),g)
        enddo
        fsigtb(i,g)=fsigt(i,g)/h_hf(i)
    enddo; enddo
endsubroutine

subroutine evaluateDtilde
    implicit none
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    
    integer i, j, k, g
    real(8)    betaL, betaR
    integer :: compL, compR
    Dtilde_c = 0.
    do g= 1, ng    
        do i=1, totm_cmfd
            ! for left boundary
            if (i.eq.1) then
                compR = i
                betaR = hDif(compR, g)/h_cmfd(i)
                Dtilde_c(1, g) = 2*albedo(1)*betaR/(albedo(1)+2*betaR)
!                Dtilde_c(1, g) = 0.
                compL = i
                compR = i+1
                betaL = hDif(compL, g)/h_cmfd(i)
                betaR = hDif(compR, g)/h_cmfd(i+1)
                Dtilde_c(2, g) = 2*betaL*betaR/(betaL+betaR)
            elseif(i.eq.totm_cmfd) then
                compL = i-1
                compR = i
                betaL = hDif(compL, g)/h_cmfd(i-1)
                betaR = hDif(compR, g)/h_cmfd(i)
                Dtilde_c(2*totm_cmfd-1, g) = 2*betaL*betaR/(betaL+betaR)
                compL = i
                betaL = hDif(compL, g)/h_cmfd(i)
                Dtilde_c(2*totm_cmfd, g) = 2*albedo(2)*betaL/(albedo(2)+2*betaL)
!                Dtilde_c(2*totm_cmfd, g) = 0.
            else
                compL = i-1
                compR = i
                betaL = hDif(compL, g)/h_cmfd(i-1)
                betaR = hDif(compR, g)/h_cmfd(i)
                Dtilde_c(2*i-1, g) = 2*betaL*betaR/(betaL+betaR)
                compL = i
                compR = i+1
                betaL = hDif(compL, g)/h_cmfd(i)
                betaR = hDif(compR, g)/h_cmfd(i+1)
                Dtilde_c(2*i , g) = 2*betaL*betaR/(betaL+betaR)
            end if
        end do
    end do
#ifdef _PRTDTLD
    write(fdtilde, '(2000(e11.5,1x))') (Dtilde_c(i*2-1,1),i=2,totm_cmfd)
#endif
end subroutine

!subroutine condense(add, wnd, wrt)
!    use define
!    use mvwnd
!    implicit none
!    include "pspec.FI"
!    include "arrays.FI"
!    include "accm.FI"
!    include "xsec.FI"
!    integer :: g, gp, i, i_f, i_c
!    real(8) :: rtot, rrmv, rsct(ng), volflux, rvolflux
!    logical,intent(in) :: add, wnd, wrt
!    
!    ! flux condense
!    phi_c = 0.
!    phi_cp = 0.
!    hTot=0;hRmv=0;hmSct=0;hNuF=0
!    
!    ! cycle flux are used to obtain initial guess
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            phi_cp(i_c, g) = phi_cp(i_c, g) + phi_f(i,g)*h(i)
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            phi_cp(i, g) = phi_cp(i, g)/h_cmfd(i)
!        end do
!    end do
!    ! surface current condense    
!    do g=1, ng
!        i_f = 1
!        do i=1, totmesh*2
!            if (jcvt2cmfd(i)) then
!                jc(i_f, g) = jf(i, g)
!                i_f = i_f + 1
!            end if
!        end do
!    end do
!    if (wrt) then
!        6(fcp, '(2000(e11.5,1x))') ((phi_cp(i,g),i=1,totm_cmfd), g=1, ng)
!        write(fcj, '(2000(e11.5,1x))') ((jc(i,g),i=1,totm_cmfd*2), g=1, ng)
!    endif 
!    
!    if (add) then
!        call add_mvwnd(phi_f, jc)
!    endif
!
!    if (wnd) then
!        call get_mvwnd(phiftemp, jc)
!    else
!       ! surface current condense    
!        do g=1, ng
!            i_f = 1
!            do i=1, totmesh*2
!                if (jcvt2cmfd(i)) then
!                    jc(i_f, g) = jfAC(i, g)
!                    i_f = i_f + 1
!                end if
!            end do
!        end do
!        phiftemp = phifAC
!    end if 
!    
!     
!    do g=1, ng
!!        rtot = 0;rrmv=0;rsct=0;
!        do i=1, totmesh
!            i_c = fi2co(i)
!            
!            ! accumulated values are used for constructing the system
!            volflux = phiftemp(i,g)*h(i)
!            phi_c(i_c, g) = phi_c(i_c, g) + volflux
!            hTot(i_c,g) = hTot(i_c,g) + xTot(complist(i),g)*volflux
!            hRmv(i_c,g) = hRmv(i_c,g) + xRmv(complist(i),g)*volflux
!            hNuF(i_c,g) = hNuF(i_c,g) + xNuf(complist(i),g)*volflux
!            do gp=1, ng
!                hmSct(g,gp,i_c) = hmSct(g,gp,i_c) + XsctM(g,gp,complist(i))*volflux
!            enddo            
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            if (phi_c(i, g)>0.) then
!                rvolflux = 1./phi_c(i, g)
!                hTot(i,g)=hTot(i,g)*rvolflux
!                hRmv(i,g)=hRmv(i,g)*rvolflux
!                hNuF(i,g)=hNuF(i,g)*rvolflux
!!#ifdef _LANISO
!!            hDif(i,g)=1./(3*(hTot(i,g)-0.1_8))            
!!#else
!!            hDif(i,g)=1./(3*hTot(i,g))
!!#endif
!                hDif(i,g)=1./(3*hTot(i,g))
!                do gp=1, ng
!                    hmSct(g,gp,i) = hmSct(g,gp,i)*rvolflux
!                end do
!                phi_c(i, g) = phi_c(i, g)/h_cmfd(i)
!            endif
!        end do
!    end do
!    
!
!    ! generate fission source
!    psi_cp = 0.
!    do g=1, ng
!        do i=1, totm_cmfd
!            psi_cp(i) = psi_cp(i) + hNuf(i, g)*phi_cp(i, g)
!        end do
!    end do        
!end subroutine


subroutine evaluateDhat
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, g

    Dhat = 0.
    do g= 1, ng    
        do i=1, totm_cmfd
            ! for left boundary
            if (hphi(i, g) .eq. 0._8) cycle
            if (i.eq.1) then 
                if (hj(i*2-1, g)==0. .and. hphi(i, g)==0) then
                    Dhat(i*2-1, g) = Dtilde_c(i*2-1, g)
                else
                    Dhat(i*2-1, g) = hj(i*2-1, g)/hphi(i, g)
                end if
                Dhat(i*2, g) = -(hj(i*2, g)+Dtilde_c(i*2, g)*(hphi(i+1, g)-hphi(i, g)))/(hphi(i+1, g)+hphi(i, g));
            elseif(i.eq.totm_cmfd) then
                Dhat(i*2-1, g) = -(hj(i*2-1, g)+Dtilde_c(i*2-1, g)*(hphi(i-1, g)-hphi(i, g)))/(hphi(i-1, g)+hphi(i, g));
                if (hj(i*2, g)==0. .and. hphi(i, g)==0) then
                    Dhat(i*2, g) = Dtilde_c(i*2, g)
                else
                   Dhat(i*2, g) = hj(i*2, g)/hphi(i, g)
                end if 
            else
                Dhat(i*2-1, g) = -(hj(i*2-1, g)+Dtilde_c(i*2-1, g)*(hphi(i-1, g)-hphi(i, g)))/(hphi(i-1, g)+hphi(i, g));
                Dhat(i*2, g) = -(hj(i*2, g)+Dtilde_c(i*2, g)*(hphi(i+1, g)-hphi(i, g)))/(hphi(i+1, g)+hphi(i, g));
            end if
        end do
    end do
    
    if (NMVWND .eq. 1) then
    do g= 1, ng    
        do i=1, totm_cmfd
            if (hphi(i, g) .eq. 0._8) cycle
            xsqsum(i,g)=xsqsum(i,g)+hphi(i,g)**2
            xsum(i,g)=xsum(i,g)+hphi(i,g)
            xysum(2*i-1,g)=xysum(2*i-1,g)+hphi(i,g)*hj(i*2-1, g)
            xysum(2*i,g)=xysum(2*i,g)+hphi(i,g)*hj(i*2, g)
            ysum(2*i-1,g)=ysum(2*i-1,g)+hj(i*2-1, g)
            ysum(2*i,g)=ysum(2*i,g)+hj(i*2, g)
        end do
    end do
    endif
!    write(fdhat, '(2000(e11.5,1x))') (Dhat(i*2,1),i=1,totm_cmfd-1)
!    write(fj, '(2000(e11.5,1x))') (hj(2*i,1),i=1,totm_cmfd-1)
!    write(fdif, '(2000(e11.5,1x))') (-(hj(2*i,1))/(hphi(i+1,1)-hphi(i,1))*h_cmfd(i),i=1,totm_cmfd-1)

!    write(fphiac, '(2000(e11.5,1x))') (hphi(i,1),i=1,totm_cmfd)
!    write(fdhat+1, '(2000(e11.5,1x))') (Dhat(i*2,1),i=1,totm_cmfd)
!    write(fdhat+2, '(2000(e11.5,1x))') (Dhat(i,1),i=1,totm_cmfd*2)
    cntacc = cntacc+1
end subroutine

#ifdef CMFDLARSEN
subroutine evaluateDhat_Larsen
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, g, ic, icc
    integer :: im1, im2
    
    logical :: jwavg=.false. ! current weighted average cross section 
    real(8) :: hjtot(totm_cmfd, ng)
    real(8) :: dtildeL2(0:totm_cmfd, 1:ng)
    real(8) :: rtot(totm_cmfd,ng), voltot(totm_cmfd,ng)
    real(8) :: chk_dif(0:totm_cmfd), chk_cur(0:totm_cmfd), chk_difrel(0:totm_cmfd)
    real(8) :: sqsumdif
    real(8) :: a(0:totm_cmfd),b(0:totm_cmfd),c(0:totm_cmfd),d(0:totm_cmfd)
!    im1=75
!    im2=88
    
!20  format    (4(es12.5,1x, 1x))
!    write(1234, 20) hj(im1*2,1), hj(im2*2,1), phi2Ls(im1), phi2Ls(im2)


    
    if (NMVWND .eq. 1) then
        phi1Lah=0
        phi2Lsh=0
        phi1Laacc=0
        htjr=0
    endif
    
    do g=1, ng
        dtildeL(0,g) = 1._8/3/htot(1,g)/delx_2(0)
        do i=1, totm_cmfd-1
            dtildeL(i,g) =(h_cmfd(i)+h_cmfd(i+1))/(htot(i,g)*h_cmfd(i) + htot(i+1,g)*h_cmfd(i+1))/3/delx_2(i) 
        enddo
        dtildeL(i,g) = 1._8/3/htot(i,g)/delx_2(i)
    enddo
    
    
!    hjtot = 0
    ! homogenization for GCMFD (additional terms)
    do g=1, ng
        phi2Lsh(0,g)=phi2Lsh(0,g)+phi2Ls(0,g)
        ic = 1
        do i=1, totmesh
            if (jcvt2cmfd(2*i)) then
                phi2Lsh(ic,g)=phi2Lsh(ic,g)+phi2Ls(i,g)
                ic = ic + 1
            end if
            icc = fi2co(i)
            phi1Laacc(icc,g)=phi1Laacc(icc,g)+phi1La(i,g)*h(i)
            htjr(icc,g)=htjr(icc,g)+xtot(complist(i),g)*phi1La(i,g)*h(i)
        end do
    end do
    do g=1,ng
        do i=1, totm_cmfd  
            hjtot(i,g) = htjr(i,g)/phi1Laacc(i,g)
            phi1Lah(i,g) = phi1Laacc(i,g)/h_cmfd(i)
        enddo
    enddo

!    if (cmfdtype .eq. 2) jwavg = .true.

    select case(cmfdtype)
    case (1)
!    if (jwavg) then ! GCMFD with Current weighted cross sections
        
        do g=1, ng
            dtildeL2(0,g) = 1._8/3/hjtot(1,g)/delx_2(0)
            do i=1, totm_cmfd-1
                dtildeL2(i,g) =(h_cmfd(i)+h_cmfd(i+1))/(hjtot(i,g)*h_cmfd(i) + hjtot(i+1,g)*h_cmfd(i+1))/3/delx_2(i) 
            enddo
            dtildeL2(i,g) = 1._8/3/hjtot(i,g)/delx_2(i)
        enddo
        
        do g=1, ng
            phi1Ls(0,g)=phi1Lah(1,g)
            do i=1, totm_cmfd-1
!                phi1Ls(i,g) =(hjtot(i,g)*h_cmfd(i)*phi1Lah(i,g) + hjtot(i+1,g)*h_cmfd(i+1)*phi1Lah(i+1,g)) \
                phi1Ls(i,g) =(htjr(i,g)+htjr(i+1,g)) \
                            /(hjtot(i,g)*h_cmfd(i) + hjtot(i+1,g)*h_cmfd(i+1)) 
            enddo 
            phi1Ls(i,g)=phi1Lah(i,g)
        enddo
    case (2)

        dtildeL2(0:totm_cmfd,1:ng)=dtildeL(0:totm_cmfd,1:ng)
        
        do g=1, ng
            phi1Ls(0,g)=phi1Lah(1,g)
            do i=1, totm_cmfd-1
!                phi1Ls(i,g) =(hjtot(i,g)*h_cmfd(i)*phi1Lah(i,g) + hjtot(i+1,g)*h_cmfd(i+1)*phi1Lah(i+1,g)) \
                phi1Ls(i,g) =(htjr(i,g)+htjr(i+1,g)) \
                            /(htot(i,g)*h_cmfd(i) + htot(i+1,g)*h_cmfd(i+1)) 
            enddo 
            phi1Ls(i,g)=phi1Lah(i,g)
        enddo
    case (3)   
        rtot=0._8;
        voltot=0._8;
        do g=1, ng
            do i=1, totmesh
                icc = fi2co(i)
                rtot(icc,g)=rtot(icc,g)+xtot(complist(i),g)*phi1La(i,g)*h(i)
                voltot(icc,g)=voltot(icc,g)+xtot(complist(i),g)*h(i)
            end do
        enddo
        
        do g=1, ng
            phi1Ls(0,g)=htjr(1,g)/voltot(1,g)
            do i=1, totm_cmfd-1
                phi1Ls(i,g)=(htjr(i,g)+htjr(i+1,g))/(voltot(i,g)+voltot(i+1,g))
            enddo
            phi1Ls(i,g)=htjr(i,g)/voltot(i,g)
        enddo
        
        do g=1, ng
            dtildeL2(0,g)= h_cmfd(1)/3/voltot(1,g)/delx_2(0)
            do i=1, totm_cmfd-1
                dtildeL2(i,g) =(h_cmfd(i)+h_cmfd(i+1))/(voltot(i,g)+voltot(i+1,g))/3/delx_2(i) 
            enddo
            dtildeL2(i,g) = h_cmfd(i)/3/voltot(i,g)/delx_2(i)
        enddo
        
        voltot=voltot/h_cmfd(1)
    end select
    

        
!    phi1ls(:,2)=0;phi2lsh(:,2)=0 ! command line for conventional CMFD
    do g=1, ng
        i=0    
        dhatl(i,g)=+dtildel(i,g)-hj(1,g)/hphi(i+1,g)
        do i=1, totm_cmfd-1
        dhatl(i,g) = (  hj(2*i,g) - phi1ls(i,g) \
                      +dtildel(i,g)*(hphi(i+1,g)-hphi(i,g))-3._8/2*dtildel2(i,g)*(phi2Lsh(i+1,g)-phi2Lsh(i-1,g)) )\
                      /(hphi(i+1,g)+hphi(i,g))
        chk_cur(i)= hj(2*i,g)
        chk_dif(i)= hj(2*i,g) - phi1ls(i,g)
        a(i)=hj(2*i,g)
        b(i)=-phi1ls(i,g) 
        c(i)=dtildel(i,g)*(hphi(i+1,g)-hphi(i,g))
        d(i)=-3._8/2*dtildel2(i,g)*(phi2Lsh(i+1,g)-phi2Lsh(i-1,g))
        
        enddo
        dhatl(i,g)=-dtildel(i,g)+hj(2*i,g)/hphi(i,g)
        write(1001, '(2000(es12.5, 1x))'), a(0:totm_cmfd)
        write(1001, '(2000(es12.5, 1x))'), -b(0:totm_cmfd)
        write(1001, '(2000(es12.5, 1x))'), c(0:totm_cmfd)
        write(1001, '(2000(es12.5, 1x))'), -d(0:totm_cmfd)
        
!        write(1001, '(2000(es12.5, 1x))'), chk_cur(0:totm_cmfd)
!        write(1001, '(2000(es12.5, 1x))'), PHI1LS(0:totm_cmfd,g)
    enddo
    chk_difrel=chk_dif/chk_cur
    
!    write(1001, '(2000(es12.5, 1x))'), chk_dif(0:totm_cmfd)
!    write(1001, '(2000(es12.5, 1x))'), chk_difrel(0:totm_cmfd)
    write(1001, '(a)'), ""
    
    sqsumdif=0._8
    do i=0, totm_cmfd
        sqsumdif = sqsumdif+chk_dif(i)**2
    enddo
!    print*, sqsumdif
    gtally(3)=sqsumdif
    tally_log = sqsumdif
end subroutine

! 2011 02 10
subroutine mt_newcmfd 
    use LU_general
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, j, g, ic, icc, gp
    integer :: im1, im2, l
    
    logical :: jwavg=.false. ! current weighted average cross section 
    real(8) :: rtot(totm_cmfd,ng), voltot(totm_cmfd,ng)
    real(8) :: sqsumdif
    real(8) :: lamda(0:totm_cmfd, ng), phi1(0:totm_cmfd, ng)
    real(8) :: lbd(totm_cmfd,ng)
    real(8) :: dd, d2, ds, avgphi1(0:totm_cmfd, ng)
!    im1=75
!    im2=88
    
!20  format    (4(es12.5,1x, 1x))
!    write(1234, 20) hj(im1*2,1), hj(im2*2,1), phi2Ls(im1), phi2Ls(im2)


    
    if (NMVWND .eq. 1) then
        phi1Lah=0
        phi2Lsh=0
        phi1Laacc=0
        htjr=0
    endif
    
    
!    hjtot = 0
    ! homogenization for GCMFD (additional terms)
    do g=1, ng
        phi2Lsh(0,g)=phi2Lsh(0,g)+phi2Ls(0,g)
        ic = 1
        do i=1, totmesh
            if (jcvt2cmfd(2*i)) then
                phi2Lsh(ic,g)=phi2Lsh(ic,g)+phi2Ls(i,g)
                ic = ic + 1
            end if
        end do
    end do

!    if (cmfdtype .eq. 2) jwavg = .true.

    voltot=0._8;
    do g=1, ng
        do i=1, totmesh
            icc = fi2co(i)
            htjr(icc,g)=htjr(icc,g)+xtot(complist(i),g)*h(i)*phi1La(i,g)
            voltot(icc,g)=voltot(icc,g)+xtot(complist(i),g)*h(i)
        end do
    enddo
    
    
    phi1(0,:)=-hj(1,:)
    do i=1, totm_cmfd
        phi1(i,:) = hj(i*2,:)
    enddo
    
    ! LAMDA (i) = lamda(i+1/2) - lamda(i-1/2)
    lamda=0
    do g=1, ng
        do i=0, totm_cmfd
            if (i.eq.0) then
                lamda(i,g)=0 !phi1(i,g)-htjr(i+1,g)/voltot(i+1,g)
!                lamda(i,g)=phi1(i,g)-htjr(i+1,g)/voltot(i+1,g)
                avgphi1(i,g)=htjr(i+1,g)/voltot(i+1,g)
            elseif (i.eq.totm_cmfd) then
                lamda(i,g)=0 !phi1(i,g)-htjr(i,g)/voltot(i,g)
!                lamda(i,g)=phi1(i,g)-htjr(i,g)/voltot(i,g)
                avgphi1(i,g)=htjr(i,g)/voltot(i,g)
            else
                lamda(i,g)=phi1(i,g)-(htjr(i,g)+htjr(i+1,g))/(voltot(i,g)+voltot(i+1,g))
                avgphi1(i,g)=(htjr(i,g)+htjr(i+1,g))/(voltot(i,g)+voltot(i+1,g))
            endif
        enddo
    enddo
    


 
    ! construct linear system
    l=0
    phi0=0
    phi0(1:totm_cmfd,:)=hphi(:,:)
    mtfull=0;
    
    do g=1, ng; do i=1, totm_cmfd
        l=l+1
        mtfull(l,l)=hRmv(i,g)*h_cmfd(i)
        lbd(l,g)=(lamda(i,g)-lamda(i-1,g))

!        ! right
!        if (i<totm_cmfd) then
!            dd=2._8/(voltot(i+1,g)+voltot(i,g))
!!            d2=-dd*(phi2Lsh(i+1,g)+phi2Lsh(i,g))/2/hphi(i+1,g)
!            d2=-dd*1/3
!            mtfull(l+1,l)=mtfull(l+1,l)+d2
!!            mtfull(l+1,l)=mtfull(l+1,l)+dd
!            lbd(l,g)=lbd(l,g)+dd*(1._8/3*hphi(i+1,g)-(phi2Lsh(i+1,g)+phi2Lsh(i,g))/2)
!            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2)/hphi(i,g)
!!            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2+lamda(i,g))/hphi(i,g)
!        else
!!            if (albedo(2) .eq. 0._8) then
!!                ds=0._8
!!            else
!!                dd=1._8/voltot(i,g)
!!                ds=dd*(-phi2Lsh(i,g)+phi2Lsh(i-1,g))/hphi(i,g)
!                ds=phi1(i,g)/phi0(i,g)
!!            endif
!        endif
!        mtfull(l,l)=mtfull(l,l)+ds
!
!        ! left
!        if (i>1) then
!            dd=2._8/(voltot(i,g)+voltot(i-1,g))
!!            d2=-dd*(phi2Lsh(i-1,g)+phi2Lsh(i-2,g))/2/hphi(i-1,g)
!            d2=-dd*1/3
!            mtfull(l-1,l)=mtfull(l-1,l)+d2
!!            mtfull(l-1,l)=mtfull(l-1,l)+dd
!            lbd(l,g)=lbd(l,g)+dd*(1._8/3*hphi(i-1,g)-(phi2Lsh(i-1,g)+phi2Lsh(i-2,g))/2)
!            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2)/hphi(i,g)
!!            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2-lamda(i-1,g))/hphi(i,g)
!        else
!!            dd=1._8/voltot(i,g)
!!            ds=dd*(phi2Lsh(i,g)-phi2Lsh(i-1,g))/hphi(i,g)
!            ds=-phi1(i-1,g)/phi0(i,g)
!        endif
!        mtfull(l,l)=mtfull(l,l)+ds
        
        if (cmfdtype .eq. 6) then
            ! right
            if (i<totm_cmfd) then
                dd=2._8/(voltot(i+1,g)+voltot(i,g))
                d2=-dd*1/3
                mtfull(l+1,l)=mtfull(l+1,l)+d2
                lbd(l,g)=lbd(l,g)+dd*(( 1._8/3-(phi2Lsh(i+1,g)+phi2Lsh(i,g))/(2*hphi(i+1,g)) )*hphi(i+1,g) \
                                     -( 1._8/3-(phi2Lsh(i,g)+phi2Lsh(i-1,g))/(2*hphi(i,g))   )*hphi(i,g)   \
                                     )
    !            lbd(l,g)=lbd(l,g)+dd*(1._8/3*hphi(i+1,g)-(phi2Lsh(i+1,g)+phi2Lsh(i,g))/2)
    !            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2)/hphi(i,g)
                ds=-d2
            else
                    ds=phi1(i,g)/phi0(i,g)
            endif
            mtfull(l,l)=mtfull(l,l)+ds
            
            ! left
            if (i>1) then
                dd=2._8/(voltot(i,g)+voltot(i-1,g))
                d2=-dd*1/3
                mtfull(l-1,l)=mtfull(l-1,l)+d2
                lbd(l,g)=lbd(l,g)-dd*(( 1._8/3-(phi2Lsh(i,g)+phi2Lsh(i-1,g))  /(2*hphi(i,g))   )*hphi(i,g)    \
                                     -( 1._8/3-(phi2Lsh(i-1,g)+phi2Lsh(i-2,g))/(2*hphi(i-1,g)) )*hphi(i-1,g)  \
                                     )
    !            lbd(l,g)=lbd(l,g)+dd*(1._8/3*hphi(i-1,g)-(phi2Lsh(i-1,g)+phi2Lsh(i-2,g))/2)
    !            ds=(dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/2)/hphi(i,g)
                ds=-d2
            else
                ds=-phi1(i-1,g)/phi0(i,g)
            endif
            mtfull(l,l)=mtfull(l,l)+ds
        else
            if (i<totm_cmfd) then
                dd=2._8/(voltot(i+1,g)+voltot(i,g))
                d2=-dd*(phi2Lsh(i+1,g)+phi2Lsh(i,g))/(2*hphi(i+1,g))
                mtfull(l+1,l)=mtfull(l+1,l)+d2
                ds=dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))/(2*hphi(i,g)) 
            else
                ds=phi1(i,g)/phi0(i,g)
            endif
            mtfull(l,l)=mtfull(l,l)+ds
            
            ! left
            if (i>1) then
                dd=2._8/(voltot(i,g)+voltot(i-1,g))
                d2=-dd*(phi2Lsh(i-1,g)+phi2Lsh(i-2,g))/(2*hphi(i-1,g))
                mtfull(l-1,l)=mtfull(l-1,l)+d2
                ds=dd*(phi2Lsh(i,g)+phi2Lsh(i-1,g))  /(2*hphi(i,g))
            else
                ds=-phi1(i-1,g)/phi0(i,g)
            endif
            mtfull(l,l)=mtfull(l,l)+ds
        endif
        

        if (cmfdtype .eq. 4) then
            mtfull(l,l)=mtfull(l,l)+lbd(l,g)/hphi(i,g)
        elseif (cmfdtype .eq. 5 .or. cmfdtype .eq. 6) then
            if ( lbd(l,g)>0) then
                mtfull(l,l)=mtfull(l,l)+lbd(l,g)/hphi(i,g)
            else
                if (l>1) then
                    if(l<totm_cmfd)then
                        mtfull(l-1,l)=mtfull(l-1,l)+lbd(l,g)/(hphi(i-1,g)+hphi(i+1,g))
                        mtfull(l+1,l)=mtfull(l+1,l)+lbd(l,g)/(hphi(i-1,g)+hphi(i+1,g))
                    else
    !                    mtfull(l-1,l)=lbd(l,g)/(2*hphi(i-1,g))
                    mtfull(l,l)=mtfull(l,l)+lbd(l,g)/hphi(i,g)
                    endif
                else
    !                mtfull(l+1,l)=lbd(l,g)/(2*hphi(i+1,g))
                    mtfull(l,l)=mtfull(l,l)+lbd(l,g)/hphi(i,g)
                endif
            endif
    !        mtfull(l,l)=hRmv(i,g)*h_cmfd(i)
        endif

        do gp=1, ng
            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
        end do
        if (i>1) mt_cmfd(1,i,g)=mtfull(l-1,l)
        mt_cmfd(2,i,g)=mtfull(l,l)
        if (i<totm_cmfd) mt_cmfd(3,i,g)=mtfull(l+1,l)
    enddo;enddo
       
    ! transpose
!    do i=1, totm_cmfd*ng
!        do j=1, totm_cmfd*ng
!            mtflu(i,j) = mtflu(j,i);
!        enddo
!    enddo
    
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)

end subroutine
#endif

!subroutine matrix_gcmfd3
!    use LU_general
!    implicit none
!    include "xsec.FI"
!    include "pspec.FI"
!    include "arrays.FI"
!    include "hmgz.FI"
!    include "CMFD_LARSEN.FI"
!    integer:: i, j, g, ic, icc, gp, tc, tcm, tcp
!    integer :: im1, im2, l
!    
!    real(8) :: fphi1(0:totm_cmfd,ng)
!    real(8) :: fsigtphi1(0:totm_cmfd,ng)
!    real(8) :: gsigtphi0(0:totm_cmfd+1,ng), gsigssphi0(0:totm_cmfd+1,ng)
!    real(8) :: lamda(0:totm_cmfd, ng), phi1(0:totm_cmfd, ng)
!    real(8) :: avgphi1(0:totm_cmfd,ng), avgtphi1(0:totm_cmfd,ng)
!    real(8) :: d, sct, l1, l2, lambda, dt, dht, gx
!    real(8) :: phi0(0:totm_cmfd+1,ng), gphi0(0:totm_cmfd+1,ng), phi2(0:totm_cmfd+1,ng)
!    real(8) :: balance(totm_cmfd,ng), source(totm_cmfd,ng), rmv(totm_cmfd,ng), leak(totm_cmfd,ng)
!    
!    ! initialize
!    fphi1=0;fsigtphi1=0;gsigtphi0=0;gsigssphi0=0;gHF=0;
!    phi0=0;gphi0=0;phi2=0;
!    ! generate the functionals according to coarse meshes
!    do g=1, ng
!    do i=1, totmesh
!        ic=fi2co(i)
!        phi0(ic,g)=phi0(ic,g)+phi0f(i,g)
!        phi2(ic,g)=phi2(ic,g)+phi2f(i,g)
!        tcm=complist(i-1)
!        tc=complist(i)
!        tcp=complist(i+1)
!        ! current term
!        fphi1(ic-1, g)=fphi1(ic-1, g)+fphi1_hf(1,i,g)
!        fphi1(ic, g)=fphi1(ic, g)+fphi1_hf(2,i,g)
!        
!        ! sig_t x phi1
!        fsigtphi1(ic-1, g)=fsigtphi1(ic-1, g)+xtot(tcm,g)*fphi1_hf(1,i,g)
!        fsigtphi1(ic, g)=fsigtphi1(ic, g)+xtot(tc,g)*fphi1_hf(2,i,g)
!        
!        ! gphi0
!        gphi0(ic-1,g)=gphi0(ic-1,g)+gphi0_hf(1,i,g)
!        gphi0(ic,g)=gphi0(ic,g)+gphi0_hf(2,i,g)
!        gphi0(ic+1,g)=gphi0(ic+1,g)+gphi0_hf(3,i,g)
!        
!        ! sig_t x phi0
!        gsigtphi0(ic-1,g)=gsigtphi0(ic-1,g)+xtot(tcm,g)*gphi0_hf(1,i,g)
!        gsigtphi0(ic,g)=gsigtphi0(ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
!        gsigtphi0(ic+1,g)=gsigtphi0(ic+1,g)+xtot(tcp,g)*gphi0_hf(3,i,g)
!
!        ! sig_ss x phi0
!        gsigssphi0(ic-1,g)=gsigssphi0(ic-1,g)+xsctm(g,g,tcm)*gphi0_hf(1,i,g)
!        gsigssphi0(ic,g)=gsigssphi0(ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
!        gsigssphi0(ic+1,g)=gsigssphi0(ic+1,g)+xsctm(g,g,tcp)*gphi0_hf(3,i,g)
!        
!        ! Equivalent X-section : nuF
!        gHF(ic-1,g)=gHF(ic-1,g)+xnuf(tcm,g)*gphi0_hf(1,i,g)
!        gHF(ic,g)=gHF(ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
!        gHF(ic+1,g)=gHF(ic+1,g)+xnuf(tcp,g)*gphi0_hf(3,i,g)
!        
!        ! Equivalent X-section : Scattering
!        do gp=1, ng
!            gHSCT(gp,g,ic-1)=gHSCT(gp,g,ic-1)+xsctm(gp,g,tcm)*gphi0_hf(1,i,g)
!            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
!            gHSCT(gp,g,ic+1)=gHSCT(gp,g,ic+1)+xsctm(gp,g,tcp)*gphi0_hf(3,i,g)
!        enddo
!    enddo
!    enddo
!    
!    ! for equivalent cross sections ( nuf & sct )
!    do g=1, ng
!        i=1
!        phi0(i,g)=phi0(i,g)/h_cmfd(i)        
!        phi2(i,g)=phi2(i,g)/h_cmfd(i)        
!        gHF(i,g)=gHF(i,g)/phi0(i,g)/h_cmfd(i)
!        do gp=1, ng
!            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)/phi0(i,g)/h_cmfd(i)
!        enddo
!
!        do i=2, totm_cmfd-1
!            phi0(i,g)=phi0(i,g)/h_cmfd(i)        
!            phi2(i,g)=phi2(i,g)/h_cmfd(i)        
!            gHF(i,g)=gHF(i,g)/phi0(i,g)/h_cmfd(i)
!            do gp=1, ng
!                gHSCT(gp,g,ic)=gHSCT(gp,g,ic)/phi0(i,g)/h_cmfd(i)
!            enddo
!        enddo
!        phi0(i,g)=phi0(i,g)/h_cmfd(i)        
!        phi2(i,g)=phi2(i,g)/h_cmfd(i)        
!        gHF(i,g)=gHF(i,g)/phi0(i,g)/h_cmfd(i)
!        do gp=1, ng
!            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)/phi0(i,g)/h_cmfd(i)
!        enddo
!    enddo
!    
!    do g=1, ng
!    do i=1, totm_cmfd-1
!        avgphi1(i,g)=(fphi1(i,g)+fphi1(i+1,g))/h_hf(i)
!        avgtphi1(i,g)=(fsigtphi1(i,g)+fsigtphi1(i+1,g))/sigth_hf(i,g)
!    enddo
!    enddo
!    
!    
!    ! balance check
!    balance=0; source=0; rmv=0; leak=0;
!    ! GCMFD formulations for inner meshes
!    do g=1, ng
!        l=totm_cmfd*(g-1)
!        do i=1, totm_cmfd
!            l=l+1
!            ! absorption
!            mtfull(l,l)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
!            rmv(i,g)=rmv(i,g)+(gsigtphi0(i,g)-gsigssphi0(i,g))
!            
!            ! left
!            d=1._8/sigth_hf(i-1,g)
!            if (i .ne. 1) mtfull(l-1,l)=-d*phi2(i-1,g)/phi0(i-1,g)
!            mtfull(l,l)=mtfull(l,l)+d*phi2(i,g)/phi0(i,g)
!            leak(i,g)=leak(i,g)+d*(phi2(i,g)-phi2(i-1,g))
!            
!            ! right
!            d=1._8/sigth_hf(i,g)
!            if (i .ne. totm_cmfd) mtfull(l+1,l)=-d*phi2(i+1,g)/phi0(i+1,g)
!            mtfull(l,l)=mtfull(l,l)+d*phi2(i,g)/phi0(i,g)
!            leak(i,g)=leak(i,g)-d*(phi2(i+1,g)-phi2(i,g))
!            
!            ! scattering
!            do gp=1, ng
!                if (gp .ne. g) then
!                    mtfull(i+(gp-1)*totm_cmfd, l) = -gHSCT(g,gp,i)
!                    source(i,g)=source(i,g)+gHSCT(g,gp,i)*phi0(i,g)
!                endif
!!                sct=hmSct(g, gp, i)*gphi0_hf(3,i-1,g)
!!                sct=sct+hmSct(g, gp, i)*gphi0_hf(2,i,g)
!!                sct=sct+hmSct(g, gp, i)*gphi0_hf(1,i+1,g)
!!                if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -sct/phi0(i,g)
!            end do
!            
!            source(i,g)=source(i,g)+1._8/keff*gHF(i,g)*phi0(i,g)
!            
!            ! lambda
!            l1=fphi1(i-1,g)-fsigtphi1(i-1,g)
!            l2=fphi1(i,g)-fsigtphi1(i,g)
!            lambda=l2-l1
!            mtfull(l,l)=mtfull(l,l)+lambda/phi0(i,g)
!
!            leak(i,g)=leak(i,g)+lambda
!            balance(i,g)=(leak(i,g)+rmv(i,g)-source(i,g))/phi0(i,g)
!        enddo
!    enddo
!    
!!    ! Typical CMFD formulation for boundary meshes
!!    do g=1, ng
!!        i=1
!!        l=totm_cmfd*(g-1)+i
!!        mtfull(l,l)=hRmv(i, g)*h_cmfd(i)
!!        ! left
!!        mtfull(l,l)=mtfull(l,l)+hj(i*2-1,g)/hphi(i,g)
!!        ! right
!!        dt=1._8/3/sigth_hf(i,g)
!!        dht=-(hj(i*2, g)+dt*(hphi(i+1, g)-hphi(i, g)))/(hphi(i+1, g)+hphi(i, g))
!!        mtfull(l+1,l)=-(dt+dht)
!!        mtfull(l,l)=mtfull(l,l)+(dt-dht)
!!        ! scattering
!!        do gp=1, ng
!!            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)        
!!        enddo
!!        
!!        i=totm_cmfd
!!        l=totm_cmfd*(g-1)+i
!!        mtfull(l,l)=hRmv(i, g)*h_cmfd(i)
!!        ! left
!!        dt=1._8/3/sigth_hf(i-1,g)
!!        dht=-(hj(i*2-1, g)+dt*(hphi(i-1, g)-hphi(i, g)))/(hphi(i-1, g)+hphi(i, g))
!!        mtfull(l,l)=mtfull(l,l)+(dt+dht)
!!        mtfull(l-1,l)=-(dt-dht)
!!        ! right
!!        mtfull(l,l)=mtfull(l,l)+hj(i*2, g)/hphi(i,g)
!!        ! scattering
!!        do gp=1, ng
!!            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)        
!!        enddo
!!    enddo
!    mtflu = mtfull
!    call LU_factorize(totm_cmfd*NG, mtflu)
!    
!end subroutine

subroutine matrix_gcmfd2
    use LU_general
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, j, g, ic, icc, gp, tc, tcm, tcp
    integer :: im1, im2, l
    
    real(8) :: fphi1(0:totm_cmfd,ng), fsigtphi1(0:totm_cmfd,ng)
    real(8) :: hc, lbd, locpsi
    real(8) :: lbdp, lbdm
    
    real(8), dimension(totm_cmfd, ng) :: balance,source,rmv,leak
    
    real(8) :: CFjh(0:totm_cmfd,ng)
    real(8) :: phi1(0:totm_cmfd, ng)
    real(8) :: dt2(totm_cmfd, ng), dh2(0:totm_cmfd, ng)
    
    real(8) :: JR, JL
    
    hc=h_cmfd(1)
    
    
    ! initialize
    fphi1=0;fsigtphi1=0

    do g=1, ng
!        phi1(0,g)=-hj(1,g)
        do i=0, totm_cmfd
            phi1(i,g)=phi1c(i,g)
        enddo
    enddo

    do g=1, ng
        do i=1, totm_cmfd-1
            dt2(i,g)= 2._8/3/hc/(htot(i,g)+htot(i+1,g))
        enddo
    enddo
    
    ! f-weighted 1st moment of total reaction =======================================================
    do g=1, ng
    do i=1, totmesh
        ic=fi2co(i)
        tcm=complist(i-1)
        tc=complist(i)
        tcp=complist(i+1)

        ! sig_t x phi1
        fsigtphi1(ic-1, g)=fsigtphi1(ic-1, g)+xtot(tc,g)*fphi1_hf(1,i,g)
        fsigtphi1(ic, g)=fsigtphi1(ic, g)+xtot(tc,g)*fphi1_hf(2,i,g)
    enddo
    enddo
    do g=1, ng
        do i=1, totm_cmfd-1
            CFjh(i,g)=(phi2(i+1,g)-phi2(i,g))/hc + (fsigtphi1(i,g))
!            CFjh(i,g)=CFjh(i,g)/(sigth_hf(i,g))
            CFjh(i,g)=CFjh(i,g)/(fsigt(i,g))
        enddo
    enddo
!    CFjh=0
    ! ===============================================================================================
    
    do g=1,ng
        i=0
        dh2(i,g)=phi1(i,g)*hc/phi0(i+1,g)
        do i=1, totm_cmfd-1
            dh2(i,g)=-(phi1(i,g)-CFjh(i,g)+dt2(i,g)*(phi0(i+1,g)-phi0(i,g))/hc)/((phi0(i+1,g)+phi0(i,g))/hc)
        enddo
        dh2(i,g)=phi1(i,g)*hc/phi0(i,g)
    enddo

    ! balance check
    balance=0; source=0; rmv=0; leak=0;
    
    ! GCMFD formulations for inner meshes
    mtfull=0._8;
    do g=1, ng
        l=totm_cmfd*(g-1)
        do i=1, totm_cmfd
            l=l+1
            mtfull(l,l) = hrmv(i,g)*hc
            ! right
            if (i .ne. totm_cmfd) then
                mtfull(l+1,l)=-dt2(i,g)-dh2(i,g)
                mtfull(l,l)=mtfull(l,l)+(dt2(i,g)-dh2(i,g))
            else
                mtfull(l,l)=mtfull(l,l)+dh2(i,g)
            endif
            
            ! left
            if (i .ne. 1) then
                mtfull(l-1,l)=-dt2(i-1,g)+dh2(i-1,g)
                mtfull(l,l)=mtfull(l,l)+(dt2(i-1,g)+dh2(i-1,g))
            else
                mtfull(l,l)=mtfull(l,l)-dh2(i-1,g)
            endif
            
            do gp=1, ng
                if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
            end do
        enddo
    enddo
    
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
    
    
    ! write balance
    write(fblc, '(2000(e12.5,1x))') (balance(i,1),i=1,totm_cmfd)
end subroutine


subroutine matrix_gcmfd3
    use LU_general
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, j, g, ic, icc, gp, tc, tcm, tcp
    integer :: im1, im2, l
    
    real(8) :: fphi1(0:totm_cmfd,ng)
    real(8) :: fsigtphi1(0:totm_cmfd,ng)
    real(8) :: gsigtphi0(0:totm_cmfd+1,ng), gsigssphi0(0:totm_cmfd+1,ng)
    real(8) :: lambda(0:totm_cmfd, ng), phi1(0:totm_cmfd, ng)
    real(8) :: avgphi1(0:totm_cmfd,ng), avgtphi1(0:totm_cmfd,ng)
    real(8) :: d, sct, l1, l2, dt, dht, gx
    real(8) :: gphi0(0:totm_cmfd+1,ng)
    real(8) :: balance(totm_cmfd,ng), source(totm_cmfd,ng), rmv(totm_cmfd,ng), leak(totm_cmfd,ng)
    real(8) :: hc, lbd, locpsi
    real(8) :: gRMV(0:totm_cmfd+1,ng)
    real(8) :: phi1a(0:totm_cmfd,ng)
    real(8) :: lbdp, lbdm
    real(8) :: gtphi0sp(3,0:totm_cmfd+1,ng), gssphi0sp(3,0:totm_cmfd+1,ng)
    real(8) :: gsctphi0sp(ng,ng,3,0:totm_cmfd+1)
    
    
    
    hc=h_cmfd(1)
    
    
    ! initialize
    fphi1=0;fsigtphi1=0;gsigtphi0=0;gsigssphi0=0;gHF=0;gHSCT=0;
    gphi0=0;
    gRMV=0;
    gtphi0sp=0;gssphi0sp=0;gsctphi0sp=0;gHF3=0;
    ! generate the functionals according to coarse meshes
    do g=1, ng
    do i=1, totmesh
        ic=fi2co(i)
        tcm=complist(i-1)
        tc=complist(i)
        tcp=complist(i+1)

        ! current term
        fphi1(ic-1, g)=fphi1(ic-1, g)+fphi1_hf(1,i,g)
        fphi1(ic, g)=fphi1(ic, g)+fphi1_hf(2,i,g)
        
        ! sig_t x phi1
        fsigtphi1(ic-1, g)=fsigtphi1(ic-1, g)+xtot(tc,g)*fphi1_hf(1,i,g)
        fsigtphi1(ic, g)=fsigtphi1(ic, g)+xtot(tc,g)*fphi1_hf(2,i,g)
        
        ! gphi0
        gphi0(ic-1,g)=gphi0(ic-1,g)+gphi0_hf(1,i,g)
        gphi0(ic,g)=gphi0(ic,g)+gphi0_hf(2,i,g)
        gphi0(ic+1,g)=gphi0(ic+1,g)+gphi0_hf(3,i,g)
        
        ! sig_t x phi0
        gsigtphi0(ic-1,g)=gsigtphi0(ic-1,g)+xtot(tc,g)*gphi0_hf(1,i,g)
        gsigtphi0(ic,g)=gsigtphi0(ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
        gsigtphi0(ic+1,g)=gsigtphi0(ic+1,g)+xtot(tc,g)*gphi0_hf(3,i,g)

        ! sig_ss x phi0
        gsigssphi0(ic-1,g)=gsigssphi0(ic-1,g)+xsctm(g,g,tc)*gphi0_hf(1,i,g)
        gsigssphi0(ic,g)=gsigssphi0(ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
        gsigssphi0(ic+1,g)=gsigssphi0(ic+1,g)+xsctm(g,g,tc)*gphi0_hf(3,i,g)

        
        ! Equivalent X-section : nuF
        gHF(ic-1,g)=gHF(ic-1,g)+xnuf(tc,g)*gphi0_hf(1,i,g)
        gHF(ic,g)=gHF(ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
        gHF(ic+1,g)=gHF(ic+1,g)+xnuf(tc,g)*gphi0_hf(3,i,g)
        
        ! Equivalent X-section : Scattering
        do gp=1, ng
            gHSCT(gp,g,ic-1)=gHSCT(gp,g,ic-1)+xsctm(gp,g,tc)*gphi0_hf(1,i,g)
            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
            gHSCT(gp,g,ic+1)=gHSCT(gp,g,ic+1)+xsctm(gp,g,tc)*gphi0_hf(3,i,g)
        enddo


        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! splitting for removal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! sig_t x phi0
        if (g_splt) then
            gtphi0sp(3,ic-1,g)=gtphi0sp(3,ic-1,g)+xtot(tc,g)*gphi0_hf(1,i,g)
            gtphi0sp(2,ic,g)=gtphi0sp(2,ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
            gtphi0sp(1,ic+1,g)=gtphi0sp(1,ic+1,g)+xtot(tc,g)*gphi0_hf(3,i,g)

            ! sig_ss x phi0
            gssphi0sp(3,ic-1,g)=gssphi0sp(3,ic-1,g)+xsctm(g,g,tc)*gphi0_hf(1,i,g)
            gssphi0sp(2,ic,g)=gssphi0sp(2,ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
            gssphi0sp(1,ic+1,g)=gssphi0sp(1,ic+1,g)+xsctm(g,g,tc)*gphi0_hf(3,i,g)
            
            ! Equivalent X-section : nuF
            gHF3(3,ic-1,g)=gHF3(3,ic-1,g)+xnuf(tc,g)*gphi0_hf(1,i,g)
            gHF3(2,ic,g)=gHF3(2,ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
            gHF3(1,ic+1,g)=gHF3(1,ic+1,g)+xnuf(tc,g)*gphi0_hf(3,i,g)
            
            ! Equivalent X-section : Scattering
            do gp=1, ng
                gsctphi0sp(gp,g,3,ic-1)=gsctphi0sp(gp,g,3,ic-1)+xsctm(gp,g,tc)*gphi0_hf(1,i,g)
                gsctphi0sp(gp,g,2,ic)=gsctphi0sp(gp,g,2,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
                gsctphi0sp(gp,g,1,ic+1)=gsctphi0sp(gp,g,1,ic+1)+xsctm(gp,g,tc)*gphi0_hf(3,i,g)
            enddo
        endif
    enddo
    enddo

    if (g_splt) then
        do g=1, ng
        do i=1, totm_cmfd
            if (i>1) gHF3(1,i,g)=gHF3(1,i,g)/phi0(i-1,g)
            if (i<totm_cmfd) gHF3(3,i,g)=gHF3(3,i,g)/phi0(i+1,g)
            gHF3(2,i,g)=gHF3(2,i,g)/phi0(i,g)
            
            do gp=1, ng
                if (i>1) gsctphi0sp(gp,g,1,i)=gsctphi0sp(gp,g,1,i)/phi0(i-1,g)
                if (i<totm_cmfd) gsctphi0sp(gp,g,3,i)=gsctphi0sp(gp,g,3,i)/phi0(i+1,g)
                gsctphi0sp(gp,g,2,i)=gsctphi0sp(gp,g,2,i)/phi0(i,g)
            enddo
        enddo
        enddo
    endif    

    
    phi1a=0;
    do g=1, ng;do i=1, totm_cmfd-1
        phi1a(i,g)=fsigtphi1(i,g)/sigth_hf(i,g)
    enddo;enddo
    
    
    ! lambda
    do g=1, ng
        i=1
        lambda(i,g)=fphi1(i,g)/hc-phi1a(i,g)
        do i=2, totm_cmfd-1
            lambda(i,g)=fphi1(i,g)/hc-phi1a(i,g)-(fphi1(i-1,g)/hc-phi1a(i-1,g))
        enddo
        lambda(i,g)=-(fphi1(i-1,g)/hc-phi1a(i-1,g))
    enddo
!    lambda=0
    
    ! for equivalent cross sections ( nuf & sct )
    do g=1, ng
        i=1
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo

        do i=2, totm_cmfd-1
            gHF(i,g)=gHF(i,g)/phi0(i,g)
            gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
            do gp=1, ng
                gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
            enddo
        enddo
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo
    enddo
    
    do g=1, ng
    do i=1, totm_cmfd-1
        avgphi1(i,g)=(fphi1(i,g)+fphi1(i+1,g))/h_hf(i)
!        avgtphi1(i,g)=(fsigtphi1(i,g)+fsigtphi1(i+1,g))/sigth_hf(i,g)
        avgtphi1(i,g)=(fsigtphi1(i,g)+fsigtphi1(i+1,g))/fsigt(i,g)
    enddo
    enddo
    
    
    ! balance check
    balance=0; source=0; rmv=0; leak=0;
    
    ! GCMFD formulations for inner meshes
    mtfull=0._8;
    do g=1, ng
        l=totm_cmfd*(g-1)
        do i=1, totm_cmfd
            l=l+1
            
            ! right
!            d=1._8/sigth_hf(i,g)
            d=1._8/fsigt(i,g)
            if (i .ne. totm_cmfd) then
                mtfull(l+1,l)=-d*phi2(i+1,g)/phi0(i+1,g) 
                mtfull(l,l)=mtfull(l,l)+d*phi2(i,g)/phi0(i,g)
                leak(i,g)=leak(i,g)-(d*phi2(i+1,g))+d*phi2(i,g)
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(2,g)/phi0(i,g)
                leak(i,g)=leak(i,g)+hc*phi1p(2,g)
            endif
            
            ! left
!            d=1._8/sigth_hf(i-1,g)
            d=1._8/fsigt(i-1,g)
            if (i .ne. 1) then
                mtfull(l-1,l)=-d*phi2(i-1,g)/phi0(i-1,g) 
                mtfull(l,l)=mtfull(l,l)+d*phi2(i,g)/phi0(i,g)
                leak(i,g)=leak(i,g)-(d*phi2(i-1,g))+d*phi2(i,g)
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(1,g)/phi0(i,g)
                leak(i,g)=leak(i,g)+hc*phi1p(1,g)
            endif
            
            if (.not.(g_splt))then
                ! absorption
                mtfull(l,l)=mtfull(l,l)+(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)*hc
            else
                mtfull(l,l)=mtfull(l,l)+(gtphi0sp(2,i,g)-gssphi0sp(2,i,g))/phi0(i,g)*hc
                if (i .ne. 1) &
                    mtfull(l-1,l)=mtfull(l-1,l)+(gtphi0sp(1,i,g)-gssphi0sp(1,i,g))/phi0(i-1,g)*hc
                if (i .ne. totm_cmfd) &
                    mtfull(l+1,l)=mtfull(l+1,l)+(gtphi0sp(3,i,g)-gssphi0sp(3,i,g))/phi0(i+1,g)*hc
            endif
            rmv(i,g)=rmv(i,g)+(gsigtphi0(i,g)-gsigssphi0(i,g))
            

            ! scattering
            do gp=1, ng
                if (.not.(g_splt))then
                    if (gp .ne. g) then
                        mtfull(i+(gp-1)*totm_cmfd, l) = mtfull(i+(gp-1)*totm_cmfd, l)-gHSCT(g,gp,i)*hc
                        source(i,g)=source(i,g)+gHSCT(g,gp,i)*phi0(i,gp)
                    endif
                else
                    if (gp .ne. g) then
                        if (i>1) mtfull(i-1+(gp-1)*totm_cmfd,l) = mtfull(i-1+(gp-1)*totm_cmfd,l)-gsctphi0sp(g,gp,1,i)*hc
                        if (i<totm_cmfd) mtfull(i+1+(gp-1)*totm_cmfd,l) = mtfull(i+1+(gp-1)*totm_cmfd,l)-gsctphi0sp(g,gp,3,i)*hc
                        mtfull(i+(gp-1)*totm_cmfd, l) = mtfull(i+(gp-1)*totm_cmfd, l)-gsctphi0sp(g,gp,2,i)*hc
                    endif
                endif
            end do
            
            locpsi=0;
            do gp=1, ng
                locpsi=locpsi+gHF(i,gp)*phi0(i,gp)
            enddo
            source(i,g)=source(i,g)+1._8/keff*locpsi*xchi(1,g)

!            ! lambda
            lbd = lambda(i,g)*hc
!            lbd =0
            lbdp=(abs(lbd)+lbd)/2
            lbdm=(abs(lbd)-lbd)/2
            mtfull(l,l)=mtfull(l,l)+lbdp/phi0(i,g)
            if (i .eq. 1) then
                mtfull(l+1,l)=mtfull(l+1,l)-lbdm/phi0(i+1,g)
            elseif(i.eq.totm_cmfd) then
                mtfull(l-1,l)=mtfull(l-1,l)-lbdm/phi0(i-1,g)
            else
                mtfull(l+1,l)=mtfull(l+1,l)-lbdm/(phi0(i+1,g)+phi0(i-1,g))
                mtfull(l-1,l)=mtfull(l-1,l)-lbdm/(phi0(i+1,g)+phi0(i-1,g))
            endif
            !mtfull(l,l)=mtfull(l,l)+lbdp/2/phi0(i,g)

!            leak(i,g)=leak(i,g)+lambda
            balance(i,g)=(leak(i,g)+rmv(i,g)-source(i,g))/phi0(i,g)
        enddo
    enddo
    
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
    
    
    ! write balance
    write(fblc, '(2000(e12.5,1x))') (balance(i,1),i=1,totm_cmfd)
end subroutine


subroutine matrix_gcmfd4 
    use LU_general
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, j, g, ic, icc, gp, tc, tcm, tcp
    integer :: im1, im2, l
    
    real(8) :: fphi1(0:totm_cmfd,ng)
    real(8) :: fsigtphi1(0:totm_cmfd,ng)
    real(8) :: gsigtphi0(0:totm_cmfd+1,ng), gsigssphi0(0:totm_cmfd+1,ng)
    real(8) :: lambda(0:totm_cmfd, ng), phi1(0:totm_cmfd, ng)
    real(8) :: avgphi1(0:totm_cmfd,ng), avgtphi1(0:totm_cmfd,ng)
    real(8) :: d, sct, l1, l2, dt, dht, gx
    real(8) :: gphi0(0:totm_cmfd+1,ng)
    real(8) :: balance(totm_cmfd,ng), source(totm_cmfd,ng), rmv(totm_cmfd,ng), leak(totm_cmfd,ng)
    real(8) :: hc, lbd, locpsi
    real(8) :: gRMV(0:totm_cmfd+1,ng)
    
    hc=h_cmfd(1)
    
    
    ! initialize
    fphi1=0;fsigtphi1=0;gsigtphi0=0;gsigssphi0=0;gHF=0;gHSCT=0;
    gphi0=0;
    gRMV=0;
    ! generate the functionals according to coarse meshes
    do g=1, ng
    do i=1, totmesh
        ic=fi2co(i)
        tcm=complist(i-1)
        tc=complist(i)
        tcp=complist(i+1)

        ! current term
        fphi1(ic-1, g)=fphi1(ic-1, g)+fphi1_hf(1,i,g)
        fphi1(ic, g)=fphi1(ic, g)+fphi1_hf(2,i,g)
        
        ! sig_t x phi1
        fsigtphi1(ic-1, g)=fsigtphi1(ic-1, g)+xtot(tc,g)*fphi1_hf(1,i,g)
        fsigtphi1(ic, g)=fsigtphi1(ic, g)+xtot(tc,g)*fphi1_hf(2,i,g)
        
        ! gphi0
        gphi0(ic-1,g)=gphi0(ic-1,g)+gphi0_hf(1,i,g)
        gphi0(ic,g)=gphi0(ic,g)+gphi0_hf(2,i,g)
        gphi0(ic+1,g)=gphi0(ic+1,g)+gphi0_hf(3,i,g)
        
        ! sig_t x phi0
        gsigtphi0(ic-1,g)=gsigtphi0(ic-1,g)+xtot(tc,g)*gphi0_hf(1,i,g)
        gsigtphi0(ic,g)=gsigtphi0(ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
        gsigtphi0(ic+1,g)=gsigtphi0(ic+1,g)+xtot(tc,g)*gphi0_hf(3,i,g)

        ! sig_ss x phi0
        gsigssphi0(ic-1,g)=gsigssphi0(ic-1,g)+xsctm(g,g,tc)*gphi0_hf(1,i,g)
        gsigssphi0(ic,g)=gsigssphi0(ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
        gsigssphi0(ic+1,g)=gsigssphi0(ic+1,g)+xsctm(g,g,tc)*gphi0_hf(3,i,g)
        
        ! Equivalent X-section : nuF
        gHF(ic-1,g)=gHF(ic-1,g)+xnuf(tc,g)*gphi0_hf(1,i,g)
        gHF(ic,g)=gHF(ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
        gHF(ic+1,g)=gHF(ic+1,g)+xnuf(tc,g)*gphi0_hf(3,i,g)
        
        ! Equivalent X-section : Scattering
        do gp=1, ng
            gHSCT(gp,g,ic-1)=gHSCT(gp,g,ic-1)+xsctm(gp,g,tc)*gphi0_hf(1,i,g)
            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
            gHSCT(gp,g,ic+1)=gHSCT(gp,g,ic+1)+xsctm(gp,g,tc)*gphi0_hf(3,i,g)
        enddo
    enddo
    enddo
    ! lambda
    do g=1, ng
        i=1
        lambda(i,g)=2._8*(1._8/hc*fphi1(i,g)-fsigtphi1(i,g)/sigth_hf(i,g) )
        do i=2, totm_cmfd-1
            lambda(i,g)=2._8*(1._8/hc*fphi1(i,g)-fsigtphi1(i,g)/sigth_hf(i,g) -(1._8/hc*fphi1(i-1,g)-fsigtphi1(i-1,g)/sigth_hf(i-1,g)) )
        enddo
        lambda(i,g)=-2._8*(1._8/hc*fphi1(i-1,g)-fsigtphi1(i-1,g)/sigth_hf(i-1,g) )
    enddo
!    lambda=0
    
    ! for equivalent cross sections ( nuf & sct )
    do g=1, ng
        i=1
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo

        do i=2, totm_cmfd-1
            gHF(i,g)=gHF(i,g)/phi0(i,g)
            gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
            do gp=1, ng
                gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
            enddo
        enddo
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo
    enddo
    
    do g=1, ng
    do i=1, totm_cmfd-1
        avgphi1(i,g)=(fphi1(i,g)+fphi1(i+1,g))/h_hf(i)
        avgtphi1(i,g)=(fsigtphi1(i,g)+fsigtphi1(i+1,g))/sigth_hf(i,g)
    enddo
    enddo
    
    
    ! balance check
    balance=0; source=0; rmv=0; leak=0;
    ! GCMFD formulations for inner meshes
    do g=1, ng
        l=totm_cmfd*(g-1)
        do i=1, totm_cmfd
            l=l+1
            ! absorption
            mtfull(l,l)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)*hc
            rmv(i,g)=rmv(i,g)+(gsigtphi0(i,g)-gsigssphi0(i,g))
            
            lbd = lambda(i,g)
            mtfull(l,l)=mtfull(l,l)+hc/2*(abs(lbd)+lbd)/phi0(i,g)
            
            ! right
            d=4._8/sigth_hf(i,g)
            if (i .ne. totm_cmfd) then
                if (i .ne. 1) then
                    mtfull(l+1,l)=-(d*phi2_hf(1,i+1,g)/phi0(i+1,g) \
                                    +hc/2*(abs(lbd)-lbd)/(phi0(i-1,g)+phi0(i+1,g)))
                else
                    mtfull(l+1,l)=-(d*phi2_hf(1,i+1,g)/phi0(i+1,g) \
                                    +hc/2*(abs(lbd)-lbd)/(phi0(i+1,g)))
                endif
                
                mtfull(l,l)=mtfull(l,l)+d*phi2_hf(2,i,g)/phi0(i,g)
                leak(i,g)=leak(i,g)-(d*phi2_hf(1,i+1,g))+d*phi2_hf(2,i,g)
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(2,g)/phi0(i,g)
                leak(i,g)=leak(i,g)+hc*phi1p(2,g)
            endif
            
            ! left
            d=4._8/sigth_hf(i-1,g)
            if (i .ne. 1) then
                if (i .ne. totm_cmfd) then
                    mtfull(l-1,l)=-(d*phi2_hf(2,i-1,g)/phi0(i-1,g) \
                                    +hc/2*(abs(lbd)-lbd)/(phi0(i-1,g)+phi0(i+1,g)))
                else
                    mtfull(l-1,l)=-(d*phi2_hf(2,i-1,g)/phi0(i-1,g) \
                                    +hc/2*(abs(lbd)-lbd)/(phi0(i-1,g)))
                endif
                mtfull(l,l)=mtfull(l,l)+d*phi2_hf(1,i,g)/phi0(i,g)
                leak(i,g)=leak(i,g)-(d*phi2_hf(2,i-1,g))+d*phi2_hf(1,i,g)
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(1,g)/phi0(i,g)
                leak(i,g)=leak(i,g)+hc*phi1p(1,g)
            endif

            ! scattering
            do gp=1, ng
                if (gp .ne. g) then
                    mtfull(i+(gp-1)*totm_cmfd, l) = -gHSCT(g,gp,i)*hc
                    source(i,g)=source(i,g)+gHSCT(g,gp,i)*phi0(i,g)
                endif
            end do
            locpsi=0;
            do gp=1, ng
                locpsi=locpsi+gHF(i,gp)*phi0(i,gp)
            enddo
!            source(i,g)=source(i,g)+1._8/keff*locpsi
            

!            leak(i,g)=leak(i,g)+lambda
            balance(i,g)=(leak(i,g)+rmv(i,g)-source(i,g))/phi0(i,g)
        enddo
    enddo
    
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
    
    ! write balance
    write(fblc, '(2000(e12.5,1x))') (balance(i,1),i=1,totm_cmfd)

end subroutine

subroutine matrix_gcmfd5
    use LU_general
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer:: i, j, g, ic, icc, gp, tc, tcm, tcp
    integer :: im1, im2, l
    
    real(8) :: fphi1(0:totm_cmfd,ng)
    real(8) :: fsigtphi1(0:totm_cmfd,ng)
    real(8) :: gsigtphi0(0:totm_cmfd+1,ng), gsigssphi0(0:totm_cmfd+1,ng)
    real(8) :: lambda(0:totm_cmfd, ng), phi1(0:totm_cmfd, ng)
    real(8) :: avgphi1(0:totm_cmfd,ng), avgtphi1(0:totm_cmfd,ng)
    real(8) :: d, sct, l1, l2, dt, dht, gx
    real(8) :: gphi0(0:totm_cmfd+1,ng)
    real(8) :: balance(totm_cmfd,ng), source(totm_cmfd,ng), rmv(totm_cmfd,ng), leak(totm_cmfd,ng)
    real(8) :: hc, lbd, locpsi
    real(8) :: gRMV(0:totm_cmfd+1,ng)
    real(8) :: phi1a(0:totm_cmfd,ng)
    real(8) :: lbdp, lbdm
    real(8) :: gtphi0sp(3,0:totm_cmfd+1,ng), gssphi0sp(3,0:totm_cmfd+1,ng)
    real(8) :: gsctphi0sp(ng,ng,3,0:totm_cmfd+1)
    
    real(8) :: dphi2(0:totm_cmfd+1,ng)
    real(8) :: tau(totm_cmfd, ng), taup(totm_cmfd,ng), taum(totm_cmfd,ng)
    real(8) :: Dj(2, 0:totm_cmfd, ng), DJTAU(2,0:totm_cmfd,ng), fphi1es(0:totm_cmfd,ng)
    real(8) :: phi1a_est(0:totm_cmfd, ng)
    
    real(8) :: mt3(3,totm_cmfd,ng)
    
    real(8) :: tau_fine(2*rcf-1, totm_cmfd, ng)
    
    integer, parameter :: methodA=1, methodS=2
    integer :: type_method
    
    if (gcmfd_opts(1) .eq. 1) then
        type_method = methodA
    else
        type_method = methodS
    endif
    
    hc=h_cmfd(1)
    
    
    ! initialize
    fphi1=0;fsigtphi1=0;gsigtphi0=0;gsigssphi0=0;gHF=0;gHSCT=0;
    gphi0=0;
    gRMV=0;
    gtphi0sp=0;gssphi0sp=0;gsctphi0sp=0;gHF3=0;
    
    dphi2 = 0;tau=0;DJ=0;taup=0;taum=0;
    ! generate the functionals according to coarse meshes
    do g=1, ng
    do i=1, totmesh
        ic=fi2co(i)
        tcm=complist(i-1)
        tc=complist(i)
        tcp=complist(i+1)

        ! current term
        fphi1(ic-1, g)=fphi1(ic-1, g)+fphi1_hf(1,i,g)
        fphi1(ic, g)=fphi1(ic, g)+fphi1_hf(2,i,g)
        
        ! sig_t x phi1
        fsigtphi1(ic-1, g)=fsigtphi1(ic-1, g)+xtot(tc,g)*fphi1_hf(1,i,g)
        fsigtphi1(ic, g)=fsigtphi1(ic, g)+xtot(tc,g)*fphi1_hf(2,i,g)
        
        ! 3D * phi2 = 1/sigt * phi2
        dphi2(ic,g)=dphi2(ic,g)+1._8/xtot(tc,g)*phi2f(i,g)
        
        
        ! gphi0
        gphi0(ic-1,g)=gphi0(ic-1,g)+gphi0_hf(1,i,g)
        gphi0(ic,g)=gphi0(ic,g)+gphi0_hf(2,i,g)
        gphi0(ic+1,g)=gphi0(ic+1,g)+gphi0_hf(3,i,g)
        
        ! sig_t x phi0
        gsigtphi0(ic-1,g)=gsigtphi0(ic-1,g)+xtot(tc,g)*gphi0_hf(1,i,g)
        gsigtphi0(ic,g)=gsigtphi0(ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
        gsigtphi0(ic+1,g)=gsigtphi0(ic+1,g)+xtot(tc,g)*gphi0_hf(3,i,g)

        ! sig_ss x phi0
        gsigssphi0(ic-1,g)=gsigssphi0(ic-1,g)+xsctm(g,g,tc)*gphi0_hf(1,i,g)
        gsigssphi0(ic,g)=gsigssphi0(ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
        gsigssphi0(ic+1,g)=gsigssphi0(ic+1,g)+xsctm(g,g,tc)*gphi0_hf(3,i,g)

        
        ! Equivalent X-section : nuF
        gHF(ic-1,g)=gHF(ic-1,g)+xnuf(tc,g)*gphi0_hf(1,i,g)
        gHF(ic,g)=gHF(ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
        gHF(ic+1,g)=gHF(ic+1,g)+xnuf(tc,g)*gphi0_hf(3,i,g)
        
        ! Equivalent X-section : Scattering
        do gp=1, ng
            gHSCT(gp,g,ic-1)=gHSCT(gp,g,ic-1)+xsctm(gp,g,tc)*gphi0_hf(1,i,g)
            gHSCT(gp,g,ic)=gHSCT(gp,g,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
            gHSCT(gp,g,ic+1)=gHSCT(gp,g,ic+1)+xsctm(gp,g,tc)*gphi0_hf(3,i,g)
        enddo


        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! splitting for removal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! sig_t x phi0
        if (g_splt) then
            gtphi0sp(3,ic-1,g)=gtphi0sp(3,ic-1,g)+xtot(tc,g)*gphi0_hf(1,i,g)
            gtphi0sp(2,ic,g)=gtphi0sp(2,ic,g)+xtot(tc,g)*gphi0_hf(2,i,g)
            gtphi0sp(1,ic+1,g)=gtphi0sp(1,ic+1,g)+xtot(tc,g)*gphi0_hf(3,i,g)

            ! sig_ss x phi0
            gssphi0sp(3,ic-1,g)=gssphi0sp(3,ic-1,g)+xsctm(g,g,tc)*gphi0_hf(1,i,g)
            gssphi0sp(2,ic,g)=gssphi0sp(2,ic,g)+xsctm(g,g,tc)*gphi0_hf(2,i,g)
            gssphi0sp(1,ic+1,g)=gssphi0sp(1,ic+1,g)+xsctm(g,g,tc)*gphi0_hf(3,i,g)
            
            ! Equivalent X-section : nuF
            gHF3(3,ic-1,g)=gHF3(3,ic-1,g)+xnuf(tc,g)*gphi0_hf(1,i,g)
            gHF3(2,ic,g)=gHF3(2,ic,g)+xnuf(tc,g)*gphi0_hf(2,i,g)
            gHF3(1,ic+1,g)=gHF3(1,ic+1,g)+xnuf(tc,g)*gphi0_hf(3,i,g)
            
            ! Equivalent X-section : Scattering
            do gp=1, ng
                gsctphi0sp(gp,g,3,ic-1)=gsctphi0sp(gp,g,3,ic-1)+xsctm(gp,g,tc)*gphi0_hf(1,i,g)
                gsctphi0sp(gp,g,2,ic)=gsctphi0sp(gp,g,2,ic)+xsctm(gp,g,tc)*gphi0_hf(2,i,g)
                gsctphi0sp(gp,g,1,ic+1)=gsctphi0sp(gp,g,1,ic+1)+xsctm(gp,g,tc)*gphi0_hf(3,i,g)
            enddo
        endif
    enddo
    enddo
    
    
    do g=1, ng
    do ic=1, totm_cmfd-1 ! surface index
        do l=1, rcf
            i=(ic-1)*rcf+l
            tc=complist(i)
            tcp=complist(i+1)
            tau_fine(l,ic, g) = tau_fine(l,ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
!            tau(ic, g) = tau(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
        enddo

        do l=0, rcf-1
            i=ic*rcf+l
            tc=complist(i)
            tcp=complist(i+1)
            tau_fine(rcf+l,ic, g) = tau_fine(rcf+l,ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
!            tau(ic, g) = tau(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
        enddo
    enddo
    enddo

    if (type_method  .eq. methodS) then
        do g=1, ng
        do ic=1, totm_cmfd-1 ! surface index
            do l=1, rcf
                i=(ic-1)*rcf+l
                tc=complist(i)
                tcp=complist(i+1)
                tau(ic, g) = tau(ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
    !            tau(ic, g) = tau(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
            enddo

            do l=0, rcf-1
                i=ic*rcf+l
                tc=complist(i)
                tcp=complist(i+1)
                tau(ic, g) = tau(ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
    !            tau(ic, g) = tau(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
            enddo
        enddo
        enddo
        do g=1, ng
        do ic=1, totm_cmfd-1
            taup(ic,g) = 1._8/2*(abs(tau(ic,g)) + tau(ic,g))
            taum(ic,g) = 1._8/2*(abs(tau(ic,g)) - tau(ic,g))
        enddo
        enddo

        
    elseif (type_method .eq. methodA) then
        do g=1, ng
        do ic=1, totm_cmfd-1 ! surface index
            do l=1, rcf-1
                i=(ic-1)*rcf+l
                tc=complist(i)
                tcp=complist(i+1)
                taup(ic, g) = taup(ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
            enddo
            do l=1, rcf-1
                i=ic*rcf+l
                tc=complist(i)
                tcp=complist(i+1)
                taum(ic, g) = taum(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
            enddo
            i=ic*rcf
            tc=complist(i)
            tcp=complist(i+1)
            if (1._8/xtot(tcp, g)>1._8/xtot(tc, g)) then
                taum(ic, g) = taum(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
                taum(ic, g) = taum(ic, g)-(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
            else
                taup(ic, g) = taup(ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(1,i,g)
                taup(ic, g) = taup(ic, g)+(1._8/xtot(tcp, g)-1._8/xtot(tc, g))*fphi2fs(2,i,g)
            endif
        enddo
        enddo
    endif
    
    
    
    do g=1, ng
        ic=1
        DJ(2,ic,g)=(dphi2(ic,g)-hc*taum(ic,g))/phi0(ic,g)             !DJ-
        do ic=2, totm_cmfd-1
            DJ(1,ic,g)=(dphi2(ic,g)-hc*taup(ic-1,g))/phi0(ic,g)       !DJ+
            DJ(2,ic,g)=(dphi2(ic,g)-hc*taum(ic,g))/phi0(ic,g)         !DJ-
        enddo
        DJ(1,ic,g)=(dphi2(ic,g)-hc*taup(ic-1,g))/phi0(ic,g)             !DJ+
    enddo
    
    do g=1, ng
        do ic=1, totm_cmfd-1
            fphi1es(ic,g) = -1_8/hc*(dj(1, ic+1,g)*phi0(ic+1,g)-dj(2, ic,g)*phi0(ic,g))
        enddo
    enddo 

    if (g_splt) then
        do g=1, ng
        do i=1, totm_cmfd
            if (i>1) gHF3(1,i,g)=gHF3(1,i,g)/phi0(i-1,g)
            if (i<totm_cmfd) gHF3(3,i,g)=gHF3(3,i,g)/phi0(i+1,g)
            gHF3(2,i,g)=gHF3(2,i,g)/phi0(i,g)
            
            do gp=1, ng
                if (i>1) gsctphi0sp(gp,g,1,i)=gsctphi0sp(gp,g,1,i)/phi0(i-1,g)
                if (i<totm_cmfd) gsctphi0sp(gp,g,3,i)=gsctphi0sp(gp,g,3,i)/phi0(i+1,g)
                gsctphi0sp(gp,g,2,i)=gsctphi0sp(gp,g,2,i)/phi0(i,g)
            enddo
        enddo
        enddo
    endif    

    
    phi1a=0;
    do g=1, ng;do i=1, totm_cmfd-1
        phi1a(i,g)=fsigtphi1(i,g)/sigth_hf(i,g)
    enddo;enddo
    
    phi1a_est=0;
    do g=1, ng
    do i=1, totm_cmfd-1
        phi1a_est(i,g)=-1._8/h_hf(i)/hc*(dphi2(i+1,g)-dphi2(i,g))
    enddo
    enddo
    
    
    
    ! for equivalent cross sections ( nuf & sct )
    do g=1, ng
        i=1
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo

        do i=2, totm_cmfd-1
            gHF(i,g)=gHF(i,g)/phi0(i,g)
            gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
            do gp=1, ng
                gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
            enddo
        enddo
        gHF(i,g)=gHF(i,g)/phi0(i,g)
        gRMV(i,g)=(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)
        do gp=1, ng
            gHSCT(gp,g,i)=gHSCT(gp,g,i)/phi0(i,g)
        enddo
    enddo
    
    do g=1, ng
    do i=1, totm_cmfd-1
        avgphi1(i,g)=(fphi1(i,g)+fphi1(i+1,g))/h_hf(i)
        avgtphi1(i,g)=(fsigtphi1(i,g)+fsigtphi1(i+1,g))/sigth_hf(i,g)
    enddo
    enddo
    
    
    ! balance check
    balance=0; source=0; rmv=0; leak=0;
    
    ! GCMFD formulations for inner meshes
    mtfull=0._8;
    do g=1, ng
        l=totm_cmfd*(g-1)
        do i=1, totm_cmfd
            l=l+1
            
            ! right
            if (i .ne. totm_cmfd) then
                mtfull(l+1,l)=-1._8/h_hf(i)*DJ(1,i+1,g) 
                mtfull(l,l)=mtfull(l,l)+1._8/h_hf(i)*DJ(2,i,g) 
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(2,g)/phi0(i,g)
            endif
            
            ! left
            if (i .ne. 1) then
                mtfull(l-1,l)=-1._8/h_hf(i-1)*DJ(2,i-1,g) 
                mtfull(l,l)=mtfull(l,l)+1._8/h_hf(i-1)*DJ(1,i,g) 
            else
                mtfull(l,l)=mtfull(l,l)+hc*phi1p(1,g)/phi0(i,g)
            endif
            
            if (.not.(g_splt))then
                ! absorption
                mtfull(l,l)=mtfull(l,l)+(gsigtphi0(i,g)-gsigssphi0(i,g))/phi0(i,g)*hc
            else
                mtfull(l,l)=mtfull(l,l)+(gtphi0sp(2,i,g)-gssphi0sp(2,i,g))/phi0(i,g)*hc
                if (i .ne. 1) &
                    mtfull(l-1,l)=mtfull(l-1,l)+(gtphi0sp(1,i,g)-gssphi0sp(1,i,g))/phi0(i-1,g)*hc
                if (i .ne. totm_cmfd) &
                    mtfull(l+1,l)=mtfull(l+1,l)+(gtphi0sp(3,i,g)-gssphi0sp(3,i,g))/phi0(i+1,g)*hc
            endif
            rmv(i,g)=rmv(i,g)+(gsigtphi0(i,g)-gsigssphi0(i,g))
            

            ! scattering
            do gp=1, ng
                if (.not.(g_splt))then
                    if (gp .ne. g) then
                        mtfull(i+(gp-1)*totm_cmfd, l) = mtfull(i+(gp-1)*totm_cmfd, l)-gHSCT(g,gp,i)*hc
                        source(i,g)=source(i,g)+gHSCT(g,gp,i)*phi0(i,gp)
                    endif
                else
                    if (gp .ne. g) then
                        if (i>1) mtfull(i-1+(gp-1)*totm_cmfd,l) = mtfull(i-1+(gp-1)*totm_cmfd,l)-gsctphi0sp(g,gp,1,i)*hc
                        if (i<totm_cmfd) mtfull(i+1+(gp-1)*totm_cmfd,l) = mtfull(i+1+(gp-1)*totm_cmfd,l)-gsctphi0sp(g,gp,3,i)*hc
                        mtfull(i+(gp-1)*totm_cmfd, l) = mtfull(i+(gp-1)*totm_cmfd, l)-gsctphi0sp(g,gp,2,i)*hc
                    endif
                endif
            end do
            
            locpsi=0;
            do gp=1, ng
                locpsi=locpsi+gHF(i,gp)*phi0(i,gp)
            enddo
            source(i,g)=source(i,g)+1._8/keff*locpsi*xchi(1,g)


            balance(i,g)=(leak(i,g)+rmv(i,g)-source(i,g))/phi0(i,g)
        enddo
    enddo
    
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
    
    
    ! write balance
    write(fblc, '(2000(e12.5,1x))') (balance(i,1),i=1,totm_cmfd)
end subroutine



subroutine matrix_shift_gcmfdPP(ke)
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"

    real(8), intent(in) :: ke
    integer :: i, j, l, m, g, gp
    real(8) :: ker
    ker = 1._8/ke
    mtflu = mtfull
    if (.not.(g_splt)) then
        do g=1, ng
            do i=1, totm_cmfd
                l = (g-1)*totm_cmfd + i
                do gp=1, ng
                    m = (gp-1)*totm_cmfd + i
                    mtflu(m,l) = mtflu(m,l)-ker*Xchi(1,g)*gHF(i,gp)*h_cmfd(i)
                end do
            end do    
        end do
    else
        do g=1, ng
            do i=1, totm_cmfd
                l = (g-1)*totm_cmfd + i
                do gp=1, ng
                    m = (gp-1)*totm_cmfd + i
                    if (i>1) mtflu(m-1,l) = mtflu(m-1,l)-ker*Xchi(1,g)*gHF3(1,i,gp)*h_cmfd(i)
                    if (i<totm_cmfd) mtflu(m+1,l) = mtflu(m+1,l)-ker*Xchi(1,g)*gHF3(3,i,gp)*h_cmfd(i)
                    mtflu(m,l) = mtflu(m,l)-ker*Xchi(1,g)*gHF3(2,i,gp)*h_cmfd(i)
                end do
            end do    
        end do
    endif
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine

subroutine makeFS_gcmfdPF
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, g
    psi =0
    if (.not.(g_splt)) then
        do g=1, ng
            do i=1, totm_cmfd
                psi(i) = psi(i) + gHF(i,g)*phi(i,g)*h_cmfd(i);
            enddo
        enddo
    else
        do g=1, ng
            do i=1, totm_cmfd
                if (i>1) psi(i) = psi(i) + gHF3(1,i,g)*phi(i-1,g)*h_cmfd(i);
                if (i<totm_cmfd) psi(i) = psi(i) + gHF3(3,i,g)*phi(i+1,g)*h_cmfd(i);
                psi(i) = psi(i) + gHF3(2,i,g)*phi(i,g)*h_cmfd(i);
            enddo
        enddo
    endif
end subroutine


subroutine make_matrix
    use LU_tdiag
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, j, k, l, g, left, right
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    do g=1, ng
        !Dhatavg = Dhatacc(:,g)/cntacc
        Dhatavg = Dhat(:,g)
        mt_cmfd(1, 1, g) = 0.
        mt_cmfd(3, 1, g) = -(Dtilde_c(2, g)+Dhatavg(2))
        mt_cmfd(2, 1, g) = -Dhatavg(1)-Dhatavg(2)+Dtilde_c(2, g)+hRmv(1, g)*h_cmfd(1)
        do i=2, totm_cmfd-1
            left = i*2-1
            right = i*2
            mt_cmfd(1, i, g) = -(Dtilde_c(left, g)+Dhatavg(left))
            mt_cmfd(3, i, g) = -(Dtilde_c(right, g)+Dhatavg(right))
            mt_cmfd(2, i, g) = -Dhatavg(left) - Dhatavg(right) + Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g)*h_cmfd(i)
        end do
        mt_cmfd(1, totm_cmfd, g) = -(Dtilde_c(totm_cmfd*2-1, g)+Dhatavg(totm_cmfd*2-1))
        mt_cmfd(3, totm_cmfd, g) = 0.
        mt_cmfd(2, totm_cmfd, g) = -Dhatavg(totm_cmfd*2-1)-Dhatavg(totm_cmfd*2)+Dtilde_c(totm_cmfd*2-1, g) + hRmv(i, g)*h_cmfd(totm_cmfd)
    end do
    
    
    do g=1, ng
        call LU_factorize(totm_cmfd, mt_cmfd(:,:,g))
    end do
end subroutine

subroutine make_matrix_LARSEN
    use LU_tdiag
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer i, j, k, l, g, left, right
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    do g=1, ng
        do i=1, totm_cmfd
            left = i-1
            right = i
            if (i>1)          mt_cmfd(1, i, g) = -(DtildeL(left,g)+DhatL(left,g))
            if (i<totm_cmfd) mt_cmfd(3, i, g) = -(DtildeL(right,g)-DhatL(right,g))
            mt_cmfd(2, i, g) = -DhatL(left,g) + DhatL(right,g) + DtildeL(left,g) + DtildeL(right,g) + hRmv(i, g)*h_cmfd(i)
        end do
    enddo
    
    
    do g=1, ng
        call LU_factorize(totm_cmfd, mt_cmfd(:,:,g))
    end do
    
end subroutine

subroutine make_matrix_withoutcurrent
    use LU_tdiag
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, j, k, l, g, left, right
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    do g=1, ng
        do i=1, totm_cmfd
            left = (i-1)*2+1
            right = i*2
            mt_cmfd(1, i, g) = -Dtilde_c(left, g)
            mt_cmfd(3, i, g) = -Dtilde_c(right,g)
            mt_cmfd(2, i, g) =  Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g) * h_cmfd(i)
        end do
        mt_cmfd(1, 1, g) = 0.
        mt_cmfd(3, totm_cmfd, g) = 0.
    end do
    do g=1, ng
        call LU_factorize(totm_cmfd, mt_cmfd(:,:,g))
    end do
end subroutine

subroutine makeFS
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, g
    psi =0
    do g=1, ng
        do i=1, totm_cmfd
            psi(i) = psi(i) + hNuF(i,g)*phi(i,g)*h_cmfd(i);
        enddo
    enddo
end subroutine

subroutine makeQ(g, k, q)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer,intent(in) :: g
    real(8),intent(in) :: k
    real(8), dimension(:), intent(out) :: q
    integer :: i,  gp
    q =0
    do i=1, totm_cmfd
        q(i) = Xchi(1,g)*psi(i)/k
        do gp = 1, ng
            if (gp.eq.g)  then 
                continue
            else
                q(i) = q(i) + hmSct(g,gp,i)*phi(i,gp)*h_cmfd(i)
            end if 
        end do
    end do
end subroutine

subroutine prolongation
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    real(8), dimension(totm_cmfd) :: ratio, psic
    real(8) :: rtnorm
    integer :: g, i, i_f, i_c
    
    srcPDF = psi/sum(psi)
    
!    psic = 0
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            psic(i) = psic(i) + psi_f(i_c)
!        enddo
!    enddo
!    do i=1, totmesh
!        i_c = fi2co(i)
!        if (psi(i_c)>0._8) then
!            adjp(i) = psi(i_c)/psic(i_c)
!        else 
!            adjp(i) = 1.
!        end if 
!    end do
    
    
end subroutine       

subroutine err_coarse(phic, err1, err2, err3) 
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    real(8), intent(out) :: err1, err2, err3
    real(8), intent(in), pointer, dimension(:, :) :: phic
    integer :: g, i
    real(8) :: rn ! ratio of normalization
    real(8) :: dif
    if (PDFREF .eq. .false.) return

    err1 = 0;err2 = 0;err3 = 0
    rn = sum(phic)/sumref
    do g=1, ng
        do i=1, totm_cmfd
            dif = abs(refc(i,g)*rn-phic(i,g))
            dif = dif /(refc(i,g)*rn)
            err1 = err1 + (dif)**2
            if (dif > err2) err2 = dif
            if (hNuF(i, g)>0 .and. dif>err3) err3 = dif
        end do
    end do
end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!         W i e l a n d t         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine make_matrix_full
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"

    integer i, j, k, l, g, left, right, gp
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    mtfull = 0.
    do g=1, ng
        !Dhatavg = Dhatacc(:,g)/cntacc
        Dhatavg = Dhat(:,g)
        i=1
        l=(g-1)*totm_cmfd + i
        mtfull(l, l) = Dhatavg(i)-Dhatavg(i+1)+Dtilde_c(i+1, g)+hRmv(i, g)*h_cmfd(i)
        do gp=1, ng
            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
        end do
        mtfull(l+1, l) = -(Dtilde_c(i+1, g)+Dhatavg(i+1))
        do i=2, totm_cmfd-1
            l=(g-1)*totm_cmfd + i
            left = i*2-1
            right = i*2
            mtfull(l-1, l) = -(Dtilde_c(left, g)+Dhatavg(left))
            mtfull(l, l) = -Dhatavg(left) - Dhatavg(right) + Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g)*h_cmfd(i)
            do gp=1, ng
                if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
            end do
            mtfull(l+1, l) = -(Dtilde_c(right, g)+Dhatavg(right))
        end do
        l=(g-1)*totm_cmfd + i
        mtfull(l-1, l) = -(Dtilde_c(i*2-1, g)+Dhatavg(i*2-1))
        mtfull(l, l) = -Dhatavg(i*2-1)+Dhatavg(i*2)+Dtilde_c(i*2-1, g) + hRmv(i, g)*h_cmfd(i)
        do gp=1, ng
            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
        end do
    end do
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine

subroutine make_matrix_full_larsen
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    
    integer i, j, k, l, g, left, right, gp
    real(8) ::offsum
    
    mtfull = 0.
    do g=1, ng
        do i=1, totm_cmfd
            l=(g-1)*totm_cmfd + i
            left = i-1
            right = i
            if (i>1) mtfull(l-1, l) = -(DtildeL(left,g)+DhatL(left,g))
            if (i<totm_cmfd) mtfull(l+1, l) = -(DtildeL(right,g)-DhatL(right,g))
            mtfull(l, l) = -DhatL(left,g) + DhatL(right,g) + DtildeL(left,g) + DtildeL(right,g) + hRmv(i, g)*h_cmfd(i) 
            do gp=1, ng
                if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
            end do
            if (i>1) mt_cmfd(1,i,g)=mtfull(l-1,l)
            mt_cmfd(2,i,g)=mtfull(l,l)
            if (i<totm_cmfd) mt_cmfd(3,i,g)=mtfull(l+1,l)
            
        end do
    end do
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine


subroutine matrix_shift(ke)
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"

    real(8), intent(in) :: ke
    integer :: i, j, l, m, g, gp
    real(8) :: ker
    ker = 1._8/ke
    mtflu = mtfull
!    do i=1, totm_cmfd*ng
!        do j=1, totm_cmfd*ng
!            mtflu(i,j) = mtfull(j,i);
!        enddo
!    enddo
    do g=1, ng
        do i=1, totm_cmfd
            l = (g-1)*totm_cmfd + i
            do gp=1, ng
                m = (gp-1)*totm_cmfd + i
                mtflu(m,l) = mtflu(m,l)-ker*Xchi(1,g)*hNuF(i,gp)*h_cmfd(i)
            end do
        end do    
    end do
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine

subroutine makeQ_wldt(k, q)
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    real(8), intent(in) :: k
    real(8), dimension(:), intent(out) :: q
    integer :: i,  g, l, gp
    q =0
    do g=1, ng
        do i=1, totm_cmfd
            l = (g-1)*totm_cmfd + i
            q(l) = Xchi(1,g)*psi(i)/k
        end do
    end do
end subroutine


subroutine fdm_coarsemesh(initk)
    use define
    use lu_tdiag
    use powermethod
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "files.FI"
    real(8), intent(out) :: initk
    real(8) :: k, err, dr
    logical :: op
    integer :: iter,i, i_c
    logical :: add, wnd
    op = FALSE
    add = .false.
    wnd = .false.
!    call condense(add, wnd, .false.)
    call evaluateDtilde
    call make_matrix_withoutcurrent
    k=1.
    call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)

    write(*, '(2(f12.5))'), k, dr
    write(hfiter, '(2(f12.5))'), k, dr
    initk = k
!    call testc_1g(k)

    if (unifsrc) then
        psi=psi*totm_cmfd/sum(psi)
        do i=1, totmesh
            i_c = fi2co(i)
            pwadj(i) = 1._8/psi(i_c)
        enddo
        do i=1, totmesh
            pwadjr(i)=1._8/pwadj(i)
        enddo
    endif
end subroutine

subroutine cmfd(k, err, iter, lmwnd, stt)
    use define
!#ifdef CMFDLARSEN
!    use lu_tdiag
!#else
!    use lu_general
!#endif
    use powermethod
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "files.FI"
    include "hmgz.FI"
    include "stsol.FI"
    include "CMFD_LARSEN.FI"
    

    real(8) :: k, err, rerr(3), dr
    logical :: op
    integer :: iter, i, j, g, i_c
    real(8) :: ratio, rnorm
    logical :: lmwnd, add
!    logical :: ifst
    integer :: nfis
    logical, intent(in) :: stt ! statistical treatment
    
    real(8) :: relfl1(totm_cmfd), relfl2(totm_cmfd-1)
    real(8) :: psit(totm_cmfd)
    
    op = FALSE

    keff = k
!    phi_f = 1.
    if (NMVWND > 0) then
        add = .true.
    else
        add = .false.
    endif
!    call condense(add, lmwnd, .true.)

#ifdef CMFDLARSEN
    if (lgcmfd) then
        if (cmfdtype .eq. 8) then
            call gcmfd4_piece
        elseif (cmfdtype .eq. 7) then
            call gcmfd3_piece
        elseif (cmfdtype .eq. 9) then
            call gcmfd2_piece
        elseif (cmfdtype .eq. 10) then
            call gcmfd5_piece
        elseif (cmfdtype>3) then
            call gcmfd_piecenew
        else
            call gcmfd_piece
        endif
    else
        call cmfd_piece
    endif
#else

    call evaluateDtilde

    !phi_cp = phi_c
    call evaluateDhat
    !call make_matrix
    call make_matrix_full
    
!    ccpsi=0.001_8

    call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
#endif 
    call prolongation
    
    write(*, '(2(2x,i5,2x,f7.5,2x,e9.3))',advance="no"), iter, k, err
    write(hfiter, '(2(2x,i5,2x,f7.5,2x,e9.3))',advance="no"), iter, k, err
    
!    endif    
    call makeFS
    nfis =0
    do i=1, totm_cmfd   
        if (psi(i)>0) nfis=nfis+1
    enddo
    rnorm = 1._8/sum(psi)*sizeprb !nfis*h_cmfd(1)
    phi=phi*rnorm
    psi=psi*rnorm/h_cmfd(1)*nfis/totm_cmfd
    
!    stmc(1:totmesh)=stmc(1:totmesh)+phi(:,1)
!    do i=1, totm_cmfd
!        phi(i,:) = phi(i,:)*rnorm/h_cmfd(i)
!    end do
    write(fcmfd, '(f10.8, 1x, 2000(es12.5,1x))') k,(psi(i),i=1,totm_cmfd)
!    do g=1, ng
!        write(fcmfdf(g), '(f10.8, 1x, 2000(es12.5,1x))') k,(phi(i,g),i=1,totm_cmfd)
!    enddo
    psit = psi
    
#ifdef _COEFF
!    write(fcp, '(2000(e11.5,1x))') ((hphi(i,g),i=1,totm_cmfd), g=1, ng)
!    write(fcj, '(2000(e11.5,1x))') ((hj(i*2,g),i=1,totm_cmfd), g=1, ng)
!    write(fcj, '(2000(e11.5,1x))') ((mtfull(i,i),i=1,totm_cmfd), g=1, ng)
    write(fcp,          '(2000(e11.5,1x))') (hphi(i,1),i=1,totm_cmfd)
    write(fcj,          '(2000(e11.5,1x))') (hj(i*2,1),i=1,totm_cmfd)
!    write(fMTCMFD,      '(2000(e11.5,1x))') (mtfull(i,i),i=1,totm_cmfd)
!    write(fMTCMFDOFF,   '(2000(e11.5,1x))') (mtfull(i-1,i),i=2,totm_cmfd)    
    i=1
    relfl1(i) =  -mtfull(i+1,i)/mtfull(i,i)
    do i=2, totm_cmfd-1
        relfl1(i) =-( mtfull(i-1,i) + mtfull(i+1,i) )/mtfull(i,i)
        relfl2(i) = -mtfull(i+1,i)/mtfull(i,i)
    end do
    relfl1(i) = -mtfull(i-1,i)/mtfull(i,i)
!    write(fRL1CMFD,   '(2000(e11.5,1x))') (relfl1(i),i=1,totm_cmfd)    
!    write(fRL2CMFD,   '(2000(e11.5,1x))') (relfl2(i),i=2,totm_cmfd-1)    
        
#endif
!    open (1005, file = "MCMFD.txt", status="unknown")
!    do j=1, totm_cmfd
!        write(1005,      '(2000(e11.5,1x))') (mtfull(i,j),i=1,totm_cmfd)
!    end do
!    close(1005)
    if (stt) then
        if (NMVWND .eq. 1) then
            stcmfdf(0)=stcmfdf(0)+k
            stsqcmfdf(0)=stsqcmfdf(0)+k**2
            do g=1, ng
                stcmfd(0,g) = stcmfd(0,g) + k
                stsqcmfd(0,g) = stsqcmfd(0,g) + k**2
            enddo
            do i=1, totm_cmfd
!                stcmfd(i) = stcmfd(i) + phi(i,1)
!                stsqcmfd(i) = stsqcmfd(i) + phi(i,1)**2
!                stcmfdf(i) = stcmfdf(i) + 1.
                stcmfdf(i) = stcmfdf(i) + psi(i)
                stsqcmfdf(i) = stsqcmfdf(i) + psi(i)**2
                do g=1, ng
                    stcmfd(i,g) = stcmfd(i,g) + phi(i,g)
                    stsqcmfd(i,g) = stsqcmfd(i,g) + phi(i,g)**2
                enddo
            enddo
!            stcmfdf(i) = stcmfdf(i) + 1._8
            if(i>totm_cmfd) then
                i=i+1
            endif
            
        else
            stcmfdf(0) = k 
            stcmfd(0,:) = k
            do i=1, totm_cmfd
                stcmfdf(i) = psi(i)
            enddo
            do g=1, ng
                do i=1, totm_cmfd
                    stcmfd(i,g)=phi(i,g)
                enddo
            enddo
        endif
    endif
    
!    write(ft,'(2000(e11.5,1x))') (stcmfdf(i),i=0,totm_cmfd)    
!    if (stt) call add_stt
    if (unifsrc) then
        psi=psi*nfis/sum(psi)
        do i=1, totmesh
            i_c=fi2co(i)
            pwadj(i) = 1._8/psi(i_c)
        enddo
        do i=1, totmesh
            pwadjr(i)=1._8/pwadj(i)
        enddo
    endif
    
contains 
subroutine cmfd_piece
    use lu_general
        call evaluateDtilde
        call evaluateDhat
        call make_matrix_full
!        cALL matrix_gcmfd2
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd_piecenew
    use lu_general
!    use lu_tdiag
!        call make_matrix_LARSEN
!        call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)
        call mt_newcmfd
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd_piece
    use lu_general
!    use lu_tdiag
        call evaluateDhat_Larsen
!        call make_matrix_LARSEN
!        call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)
        call make_matrix_full_larsen
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd2_piece
    use lu_general
        call matrix_gcmfd2
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd3_piece
    use lu_general
!    use lu_tdiag
!        call make_matrix_LARSEN
!        call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)
        call matrix_gcmfd3
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS_gcmfdPF, makeQ_wldt, matrix_shift_gcmfdPP, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd5_piece
    use lu_general
!    use lu_tdiag
!        call make_matrix_LARSEN
!        call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)
        call matrix_gcmfd5
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS_gcmfdPF, makeQ_wldt, matrix_shift_gcmfdPP, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
subroutine gcmfd4_piece
    use lu_general
        call matrix_gcmfd4
        call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS_gcmfdPF, makeQ_wldt, matrix_shift_gcmfdPP, solve_LU, ccpsi, k, err, dr, iter)
end subroutine
endsubroutine


subroutine add_stt
    implicit none
    integer :: i
!    stcnt=stcnt+1
!    do i=1, ndim
!        sumpsi(i)=sumpsi(i)+psi(i)
!        sqsumpsi(i)=sqsumpsi(i)+psi(i)**2
!    enddo
end subroutine

subroutine print_stt
    implicit none
    integer :: i
    real(8), pointer, dimension(:) :: avg, std
!    call dmalloc(avg, ndim)
!    call dmalloc(std, ndim)
!    
!    if (stcnt>0) then
!        do i=1, ndim
!            avg(i)=sumpsi(i)/stcnt
!            std(i)=sqrt((sqsumpsi(i)-sumpsi(i)**2/stcnt)/stcnt/(stcnt-1))
!        enddo
!!        open (4001, file='avgcmfd.txt", status='unknown')
!    else
!    endif
    
end subroutine

function SHN_CMFD
    implicit none
    include "pspec.FI"
    real(8) :: SHN_CMFD, sumpsi, pj
    integer :: i, g, fismesh
	SHN_CMFD = 0.
	fismesh =0
    sumpsi = sum(psi)
!	do i=1, totmesh
!	    if (psi_f(i)>0.) fismesh = fismesh+1
!	end do
	
	do i=1, totm_cmfd
	    pj = psi(i)/sumpsi
	    if (pj>0.) SHN_CMFD = SHN_CMFD -pj*log(pj)
	end do
	SHN_CMFD = SHN_CMFD/ log(2.)
end function


subroutine testc_1g(k)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    real(8), intent(in) :: k
    integer :: i,j, o, m
    real(8) :: leak(3), absr(3), src(3), val(3), tot(3)
    real(8) :: diff(300), diff2(150), phis, phim, phio
    ! 1 group test
    call makeCurrent
    do i=1, totm_cmfd
        leak(1) = hj(i*2-1,1)+hj(i*2,1)
        leak(2) = hj(i*2-1,2)+hj(i*2,2)
        leak(3) = leak(1) + leak(2)
        absr(1) = phi(i,1)*hRmv(i, 1)*h_cmfd(i)
        absr(2) = phi(i,2)*hRmv(i, 2)*h_cmfd(i)
        src (1) = 1._8/k*psi_c(i)
        src (2) = hmSct(2,1,i)*phi(i,1)*h_cmfd(i)
        src (3) = src(1)
        absr(3) = absr(1)+absr(2)-src(2)
        tot (1) = hTot(i,1)*phi(i,1)*h_cmfd(i)
        tot (2) = hTot(i,2)*phi(i,2)*h_cmfd(i)
        tot (3) = (tot(1)+tot(2))/(phi_f(i,1)+phi(i,2))/h_cmfd(i)
        val (1) = (leak(1) + absr(1) - src(1))/phi(i,1)
        val (2) = (leak(2) + absr(2) - src(2))/phi(i,2)
        val (3) = (leak(3) + absr(3) - src(3))/(phi(i,1)+phi(i,2))
        diff2(i) = 1._8/3/tot(3)
!        if (i<150) diff(i) = -(jf(i*2,1)+jf(i*2,2))/(phi_f(i+1,1)+phi_f(i+1,2)-(phi_f(i,1)+phi_f(i,2)))*h(i)
        j = i*2-1
        if (j>1 .and. j<2*totm_cmfd) then 
            m = i
            o = i-1
            phio = phi(o,1)+phi(o,2)
            phim = phi(m,1)+phi(m,2)
            phis = (sum(hDif(o,:))*phio + sum(hdif(m,:))*phim)/(sum(hDif(o,:))+sum(hDif(m,:)))
            diff(j) = -h_cmfd(i)/2*(jf(j,1) + jf(j,2))/(phis-phim)
        endif
    end do
    contains
    
subroutine makeCurrent
    use define
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer :: g, i, i_f, i_c
    
    ! surface current condense    
    do g=1, ng
        do i=1, totm_cmfd
            if (i.eq.1) then
!            hj(i*2-1, g) = Dtilde_c(i*2-1, g) * hphi(i, g)
            hj(i*2-1, g) = Dtilde_c(1, g) * hphi(i, g)
            hj(i*2, g) = -Dtilde_c(i*2, g) * (hphi(i+1, g) - hphi(i, g))
            elseif(i.eq.totm_cmfd) then
            hj(i*2-1, g) = -Dtilde_c(i*2-1, g) * (hphi(i-1, g) - hphi(i, g))
!            hj(i*2, g) = Dtilde_c(i*2, g) * hphi(i, g)
            hj(i*2, g) = Dtilde_c(1, g) * hphi(i, g)
            else
            hj(i*2-1, g) = -Dtilde_c(i*2-1, g) * (hphi(i-1, g) - hphi(i, g))
            hj(i*2, g) = -Dtilde_c(i*2, g) * (hphi(i+1, g) - hphi(i, g))
            end if 
        end do
    end do
end subroutine
    
end subroutine 

end module
